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Abstract: We use a systematic effective field theory setup to derive the bbH production 
cross section. Our result combines the merits of both fixed 4-flavor and 5-flavor schemes. 
It contains the full 4-flavor result, including the exact dependence on the 6-quark mass, 
and improves it with a resummation of collinear logarithms of mh/rriH- In the massless 
limit, it corresponds to a reorganized 5-flavor result. While we focus on bbH production, 
our method applies to generic heavy-quark initiated processes at hadron colliders. Our 
setup resembles the variable flavor number schemes known from heavy-flavor production 
in deep-inelastic scattering, but also differs in some key aspects. Most importantly, the 
effective 6-quark PDF appears as part of the perturbative expansion of the hnal result 
where it effectively counts as an 0{as) object. The transition between the fixed-order (4- 
flavor) and resummation (5-flavor) regimes is governed by the low matching scale at which 
the 6-quark is integrated out. Varying this scale provides a systematic way to assess the 
perturbative uncertainties associated with the resummation and matching procedure and 
reduces by going to higher orders. We discuss the practical implementation and present 
numerical results for the bbH production cross section at NLO-I-NLL. We also provide a 
comparison to the corresponding predictions in the hxed 4-flavor and 5-flavor results and 
the Santander matching prescription. Compared to the latter, we find a slightly reduced 
uncertainty and a larger central value, with its central value lying at the lower edge of our 
uncertainty band. 
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1 Introduction 


The formulation of reliable predictions for heavy-quark initiated processes has been the 
subject of much study over many years, in particular for the determination of parton 
distribution functions (PDFs) in the context of deep-inelastic scattering (DIS) [1-17]. At 
the LHC, important examples of heavy-quark initiated processes are Higgs or vector-boson 
production in association with heavy quarks. In this paper, we are interested in Higgs 
production in association with b quarks, i.e., the inclusive 66F7-induced cross section. 

In a typical hard-scattering process with protons in the initial state, there are at 
least two parametrically separate scales. First, the hard scale fin ~ Q, where Q denotes 
the physical quantity that determines the momentum transfer in the hard interaction, e.g., 
in DIS or Q = mn for the case of Higgs production that we will be interested in. 
Second, the low scale /ta ~ ^.qcd, which separates the perturbative and nonperturbative 
regimes and is typically taken to be of order the proton mass, /ta ~ 1 GeV. In the limit 
/^A Q we can apply the standard QCD factorization theorem [18-20] to compute the 
hadronic cross section in terms of the partonic cross section convolved with PDFs. 

For heavy-quark initiated processes, the mass m of the heavy quark introduces another 
physical scale. Depending on its value, we can distinguish two parametrically different 
cases, shown in figure id 

(a) m ~ Q: There is a single parametric scale ~ m ~ Q in addition to 

(b) m Q: There are two parametric scales fin ~ Q and Hm ~ in addition to ^a- 

When working in the limit m ^ Q, the heavy quark never appears in the initial state of 
the hard partonic process. Instead, it is produced as part of the hard interaction at fj,H by 
an incoming gluon splitting into a pair of heavy quarks. The partonic calculation contains 
the exact dependence on m, including the correct m-dependent phase space. The gluon 
splitting into a heavy-quark pair contains a collinear singularity, which is regulated by m, 
and as a result produces logarithms ln{m/Q). For m ^ Q, these collinear logarithms are 
counted as small and are included at fixed order in the as expansion. 

When working in the limit m Q, the heavy quark explicitly appears in the initial 
state of the hard partonic process, and the collinear logarithms are resummed to all orders 
in Us into an effective heavy-quark PDF. The quark mass m only appears in the boundary 
condition of the PDF’s DGLAP evolution, which starts at the scale ~ The hard 
process itself is computed in the m —)> 0 limit. That is, hnite-mass effects of 0{jnlQ\ 
including the exact phase space of the gluon splitting into a massive quark pair, are power 
corrections and are neglected. 

Predictions obtained in strictly one of the above two limits are usually referred to as 
obtained in a fixed-flavor number scheme. Which of these limits is more appropriate in 

^In principle, there is a third parametric limit m ^ Q, which we are not interested in. In this case, 
when the heavy quark appears as an external state, m itself is the physical quantity that sets the hard 
interaction scale, so Q = m. The relevant setup is then determined by what other parametrically smaller 
physical scales are present in the process. Otherwise, when the heavy quark only appears in internal loops, 
it can simply be integrated out. 
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Figure 1. The two parametric scale hierarchies for inclusive cross sections for heavy-quark initiated 
processes. 


practice depends on the process and the numerical size of the corrections. For 6-initiated 
processes at hadron colliders, the relative importance of ln{mb/Q) corrections has been 
discussed for example in [21-23]. 

To obtain the best possible theoretical predictions, it is often desirable to have a 
complete description that incorporates the results from both limits. In this way, the final 
result is valid in each limit as well as in the transition region in between, and hence one 
can be agnostic about which parametric regime is the more appropriate one. 

For bbH, predictions exist in the 4-flavor scheme (4FS) [24, 25], which works in the 
limit mb ~ mjy, and in the 5-flavor scheme (5FS) [26-29], which works in the limit mb 
mn- Currently, both predictions are combined using the pragmatic “Santander Matching” 
prescription [30], which is a weighted average of the 4FS and 5FS predictions, where the 
relative weighting depends on the numerical size of ln(mft/mjy). 

There are various methods available in the literature, referred to as variable-flavor 
number schemes (VFNS), which aim to combine the virtues of both limits in a more 
systematic fashion. That is, they include the full m dependence in the limit m ^ Q and 
the resummation of collinear logarithms \a{m/Q) in the limit m Q. There are a number 
of such schemes available, namely the ACOT scheme [1, 2] (and its simplified variants 
S-ACOT [6], S-ACOT-x [7, 13], and the more recent m-ACOT [31] for hadron-hadron 
collisions), the TR scheme [3, 8], and the FONLL scheme [11, 32]. The differences between 
the schemes essentially amount to how the two limits are combined. 

Effective field theories (EFTs) are the standard tool to describe processes with para¬ 
metrically separated scales, allowing to systematically resum the logarithms of ratios of 
these scales. In this paper, we discuss the EFT formulation of heavy-quark initiated pro¬ 
cesses for the case of inclusive cross sections. All the basic ingredients are actually well 
known in this case. Nevertheless, we find it worthwhile to discuss the EFT formulation 
in detail, as it provides a conceptually clear field-theoretic derivation, including the tran¬ 
sition between the two parametric regimes and a way to assess the associated theoretical 
uncertainties. This setup can also be extended to more differential cross sections, which 
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we leave for future work. A similar setup has also been used to incorporate quark-mass 
effects for final-state jets in Refs. [33-35]. 

Our final result for DIS resembles the aforementioned schemes in several ways, but 
also differs in some key aspects. Most importantly, the 6-quark PDF is not treated as an 
external 0(1) quantity. Rather, it contributes as part of the perturbative series of the 
final result, where it effectively counts as an 0(as) object. (In this work, we follow the 
assumption made in all available PDF sets that there is no intrinsic bottom in the proton, 
such that the effective bottom-quark PDF is generated purely perturbatively.) 

The application of our method to hadron-hadron collisions is completely straightfor¬ 
ward. Our final result for bbH encompasses the merits of both 4F and 5F schemes. It 
contains the full 4FS result at NLO, including the exact mb dependence and phase space. 
In addition, it improves the 4FS result with the all-order resummation of collinear loga¬ 
rithms up to NLL order. In the —>■ 0 limit, our result corresponds to a reorganized 5FS 
result, where the perturbative series is expanded to NLO with the 6-quark PDF counted 
as 0{as). 

In the next section, we discuss the general setup in detail, focussing on DIS to be 
specific. In section 3 we briefly discuss the similarities and differences with respect to 
other heavy-flavor schemes in the literature. Then in section 4, we apply this framework to 
bbH production. We discuss in detail the perturbative uncertainties and present our final 
numerical results at LO-I-LL and NLO-I-NLL. We also compare to the predictions in the 
4F and 5F schemes using a consistent set of inputs. We conclude in section 5. 

2 EFT formulation of heavy-quark initiated processes 

In this section, we discuss the EFT formulation in detail. For simplicity and to be specific 
we frame the discussion in the context of heavy-quark production in DIS, where we have 
to deal with only one strongly-interacting initial state. In this case we associate Q = 
■sj —cp' ■ The extension to hadron-hadron collisions is straightforward and will be discussed 
in section 4. For definiteness we consider the heavy quark to be the 6 quark, ^ and treat 
the four lighter quarks as massless. We take Q < mt, so we can essentially ignore the top 
quark (i.e., we either integrate it out at the scale fin ~ Q or it has already been integrated 
out at a higher scale). 

In section 2.1, we review the case fi\ <C mb ~ Q, corresponding to figure 1(a), where 
a single matching step at the hard scale hh ~ Q ~ mb is required. We will refer to this 
as the fixed-order region or limit. This also serves to introduce our notation and language. 
In section 2.2, we discuss the case ha ^ mb <C Q, corresponding to figure 1(b), where two 
separate matching steps, at hh ~ Q and fim ~ mb, are performed. We will refer to this as 
the resummation region or limit. In section 2.3, we discuss the appropriate perturbative 

^Our setup can be equally applied to processes involving the top qnark. For the charm quark, the low 
value of its mass might not jnstify the treatment rric ^ Ha, which wonld mean that pa/uIc corrections are 
important. In this case, a better treatment would be to take rric ~ /ta and not integrate out the charm 
quark, but instead consider a nonperturbative charm PDF. 
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counting in our result, and in section 2.4, we combine the results in both limits to yield 
our final predictions valid in both limits and anywhere in between. 

Throughout this paper, we use roman indices i,j, fe to denote the light flavors, i.e., 
the four light quarks and the gluon. We also use the convention that any repeated indices 
are implicitly summed over (also a repeated index b implies a sum over b and b). For 
clarity, we will focus on the dependence on the relevant physical and renormalization scales, 
but suppress all other kinematic dependences. In particular, we will not write out the 
dependence on the momentum fractions and the Mellin-type convolutions in them. We 
will denote the number nj of light active flavors as superscripts for quantities where the 
distinction is relevant, e.g., vs. 

2.1 mb ~ Q: Fixed order 

In this case, shown in figure 1(a), the 6-quark mass is treated parametrically as of the 
same size as Q. At the scale fin r\j Q ~ mb, all degrees of freedom with virtualities 
~ m^, including the heavy 6 quark, are integrated out. We match full QCD onto a 
theory of collinear gluons and collinear light quarks with typical virtuality This 

matching step is precisely equivalent to the standard operator product expansion (OPE) 
in DIS [36-45], which we briefly review now. 

We dehne the DIS operator Odis(Q) 'fnb), whose proton matrix element determines the 
DIS cross section (or equivalently the hadronic tensor or DIS structure functions), 

da{Q,mb) = {p\OBis{Q,mb)\p) 

At the scale pn, h is matched onto a sum of nonlocal PDF operators 


where a sum over light quarks and gluons i = u,d,s,c,g is understood, and “(g)” de¬ 
notes the Mellin-type convolutions in the momentum fractions. The Wilson coefficients 
Di{Q,mb, ph) are also called coefficient functions. The o'f\p) are the standard MS- 
renormalized quark and gluon PDF operators [45]^, whose proton matrix elements dehne 
the nonperturbative PDFs, 

ff\p) = {p\0f\p)\p). (2.3) 

Since the 6 quark is being integrated out and not present in the theory below pn, there is 
also no Ob operator and no Db coefficient on the right-hand side of eq. (2.2). As indicated, 
the right-hand side of eq. (2.2) is the leading term in an expansion in Aq(-,j^/Q^, where the 

®In SCET, this is the purely n-collinear quark and gluon sector, which is equivalent to a boosted version 
of QCD, where n'^ = (l,n) and n is the direction of the incoming proton. In lightcone coordinates, 
the momentum of the collinear modes scales as pc ~ (Q, Aqcd/Q, hqcD). In principle, there could also 
be soft modes with momentum scaling ps ~ (Aqcd, Aqcd, Aqcd), and also Glauber modes. Since their 
contributions cancel in the inclusive cross section [20], they are not needed here. 

^For corresponding operator definitions in SCET and a discussion of their equivalence see e.g. refs. [46- 
48]. 


( 2 . 2 ) 


OBisiQ,'mb) = Di{Q,mb,pH) 




1 + 0 


^QCD 


( 2 . 1 ) 
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low scale is set by the external proton state we are eventually interested in. For ease 

of notation, we will not indicate these power corrections in the rest of this section. 

In the above, as(^) = and are renormalized with nj = 4 active quark 

flavors. That is, we use MS with dimensional regularization with respect to the four light 
quark flavors, while 6-quark loops are renormalized in the decoupling scheme, such that 
the 6-quark decouples from the theory below hh (see Appendix A). 

Since Odis determines the full-theory cross section, it does not have an explicit depen¬ 
dence on ///f, i.e., it does not receive additional operator renormalization. It only has an 
implicit dependence on through the renormalization of asifJ-n), which cancels order by 
order in perturbation theory. On the other hand, the coefficients Di are explicitly ^ depen¬ 
dent, and their /r dependence cancels against the explicit /r dependence of the operators 

The full dependence on the physical scales Q and mh, which are treated as hard scales, 
resides in the Wilson coefficients Di(Q, m;,, ^ih)- The coefficients Di at some scale /x contain 
logarithms ln(/x/Q) ~ ln(^/m6). Therefore, they are computed by a perturbative matching 
calculation (see below) at the hard scale ~ Q ~ m;,, where they contain no large 
logarithms. 

The PDFs fi{ix) at some scale contain logarithms ln(/x/AQCD)- Hence, the in¬ 
put PDFs that are determined from the experimental data are defined at a low scale 
/^A ^ AqcD) which should still be large enough for perturbation theory to be valid. All 
contributions from lower scales, including the nonperturbative regime, are absorbed into 
the input PDFs fi{nA). The renormalization of the PDF operators leads to their renor¬ 
malization group equation (RGE) 

' (6^) = Ty ® Of^ (a) , (2.4) 

where are the PDF anomalous dimensions, which are given in terms of the standard 
Altarelli-Parisi splitting functions [38, 40, 43]. The solution of this RGE yields the standard 
DGLAP evolution [43, 49, 50] relating the operators (and PDEs) at a scale /Uq to the 
operators (and PDEs) at the scale jj,, 

= (2.5) 

As denoted, the anomalous dimensions and evolution factors involve Uf = A light quark 
flavors. By definition, the coefficients and operators in eq. (2.2) must be evaluated at the 
same scale, which means the operators on the right-hand side give PDFs at containing 
large logarithms These logarithms are resummed by using eq. (2.5) to evolve 

the PDFs from the low scale up to 

/f' {^^H) = {fiH, Ma) ® /f' M ■ (2.6) 

Equivalently, we can perform the resummation for the Wilson coefficients. The coef¬ 
ficients and operators obey inverse RGEs, since their scale dependences must cancel each 
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Figure 2. Schematic leading-order matching for heavy-quark production in DIS for mj, ^ Q. 



Figure 3. Schematic NLO matching for heavy-quark production in DIS for mt ~ Q. 


other. After performing the matching at the scale hh, the coefficients are evolved from /ijy 
down to /TAj 

Dj{Q,mb,fiA) = ® ulf {^h, fiA) ■ (2.7) 

The evolution factor is precisely the same as in eq. (2.6). Evolving the coefficients down 
corresponds to successively integrating out virtualities between and After evolving 
down to fiA, we can take the proton matrix element to obtain the final DIS cross section 

da^^{Q,mb) = Di{Q,mb, hh) ^ ulf {hh, fiA) ® ■ ( 2 . 8 ) 

The full cross section on the left-hand side contains large logarithms ln((5/AQCD) 
ln(m;,/AQCD)) which on the right-hand side are factorized into logarithms \n{Q/fj,H) and 
ln(mb /which are considered small and reside in the coefficients, large logarithms 
fj,H/ a), which are resummed into the evolution factor, and logarithms ln(/XA/AQCD), 
which are absorbed into the PDFs. This result for the cross section is precisely the 4FS 
result. Since the mb dependence is included at fixed order in eq. (2.8), we will refer to it 
as the fixed-order (“FO”) result. 

2.1.1 Matching at ^ Q ^ mb 

The Wilson coefficients Di(Q,mb, fJ-n) are determined in perturbation theory by matching 
the matrix elements of both sides of eq. (2.2) between the same partonic external states j, 

{j\0-Dis{Q,'rnb)\j) = Di{Q,mb, ^j.h) ® {j\0^f^{fLH)\j) ■ (2.9) 

The left-hand side corresponds to the full-theory matrix element. The partonic matrix 
element of the MS-renormalized PDF operators on the right-hand side are the partonic 
PDFs, 

flfjM = {j\of\fiH)\j) = . (2.10) 
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They are equivalent to the collinear MS subtractions often denoted as 

The matching in eq. (2.9) is performed order by order in as{^iH) = for which 

we expand each of the pieces as 


{j\0-Dis{Q,mb)\j) = ^(j|ODis(Q,?^6)|j 




dvr 


1 ^ 


/i/j(M«)=E^/i 


[4](fc) 




1 k 


Att 




asifJ-n) 


1 k 


dvr 


( 2 . 11 ) 


The leading-order matching is shown schematically in figure 2. At the lowest order, only 
the gluon external state contributes. Light quarks in the external state first contribute 
at NLO. The 5-quark does not appear as external state in the matching calculation, since 
it cannot appear anymore on the right-hand side. Writing out the dependence on the 
momentum fraction z explicitly, the partonic PDFs at LO are simply given by 

= ( 2 . 12 ) 

so the LO gluon coefficient is directly given by the LO diagram on the left of figure 2, 

= {9\ODis{Q,mb)\g)^^\ 

Djj^\Q,mb,fiH) = 0. (2.13) 


The schematic matching at NLO is shown in figure 3. There are virtual and real 
emission corrections to the gluon channel as well as the contribution with light (anti) quarks 
in the external state. Using pure dimensional regularization to regulate both the UV and 
IR, the MS-renormalized partonic PDFs at NLO are 




^(^) = -- e{z)Pggiz) + MA) 5(1 - z) 


€ L 


f^^}^"\^) = -l‘^CF9{z)Pgg{z)^ 


(2.14) 


where /3o{nf) = (IIC^—drpn/)/3 (with n/ = d here) and the one-loop (LO) gluon splitting 
functions are 


^ggi^) — 2 


(l-z + z^) 


2\2 r 


9(1-z) 


1 - z 


Pggiz) = 9(1 - Z) 


1 + (1-Z) 


(2.15) 


Together with the LO coefficients from eq. (2.13), the NLO matching coefficients are 


Df\Q,mb,9H) = {g\ODis(Q,mb)\g)^'^^ - D^^\Q,mb, fin) ^ 

Df\Q,mb,gH) = (g|C>Dis(Q, - D^g^Q^mb, ® ■ (2.16) 
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The 1/e terms in the in eq. (2.14) are collinear IR divergences. They precisely cancel 
between the two terms on the right-hand side, such that NLO Wilson coefficients are free 
from IR divergences. While the same IR regulator must be used in the full and effective 
theories, the Wilson coefficients are independent of the specific IR regulator. On the other 
hand, the coefficients do explicitly depend on the UV renormalization scheme of the PDF 
operators, which is where the standard MS scheme is used. In other words, the fact 
that eq. (2.14) are a pure pole contribution is just an artifact from using pure dimensional 
regularization for both UV and IR divergences.^ Using any other IR regulator, e.g., putting 
the external states off shell, the would look different, but the final results for the 
would be exactly the same once the IR regulator is taken to zero. 

2.2 rrib <C Q: Resummation 

In this case, shown in figure 1(b), there is a parametric hierarchy between the 6-quark 
mass mb and Q. The calculation now proceeds via a two-step matching. First, at the scale 
~ Q, all degrees of freedom with virtualities ~ are integrated out and full QCD is 
matched onto a theory of collinear gluons, collinear light quarks, and in addition collinear 
massive 6-quarks, all with typical virtualities ~ Next, we evolve from pn down to 
the intermediate scale pm ~ fnb- At pm, all degrees of freedom with virtualities ~ are 
integrated out, including the massive 6-quark, and the theory is matched onto the same 
theory as in the previous section 2.1 of collinear gluons and collinear light quarks with 
typical virtuality 

The matching at the scale pn proceeds as before, except that above pm the bottom 
quark is still a dynamical degree of freedom. Analogous to eq. (2.2), the DIS operator 
is matched onto a sum of nonlocal PDF operators, which now includes the bottom-PDF 
operator Ob, 


Here, we use the notation Ci^b for the Wilson coefficients to distinguish them from the Di 

f5l f5l 

coefficients in the previous subsection. The PDF operators, 0\ and , have the same 
structure as those in eqs. (2.2) and (2.3). The essential difference is that they are now 
renormalized with nj = 5 active flavors. That is, 6-quark loops are now renormalized using 
MS with dimensional regularization. 

Eq. (2.17) corresponds again to the standard OPE in DIS. It is important to note, 
however, that the expansion performed in eq. (2.17) is by construction an expansion in 
where is the typical virtuality of the external states in the theory below pH. 

^Technically, all loop corrections to the bare PDF matrix elements are scaleless and vanish, which means 
the UV and IR divergences are precisely equal with opposite sign and cancel each other. Adding the UV 
counterterms then leaves the IR divergences. In this case, the /i dependence in is purely through 

oisifJ.) and the coefficients have no explicit p dependence as written in eq. (2.11). 

®In SCET this would be a theory containing massive collinear fermions [51]. The collinear modes 
have momentum scaling pc ~ {Q,ml/Q,mb). The corresponding soft modes with momentum scaling 
Ps ~ {mb, mb, mb) are again not needed since they cancel. This implies that the production of secondary 
6-quarks can only arise from the splitting of collinear gluons. 


(2.17) 


OBis{Q,mb) = Ci{Q,pH) ® of\pH) + Cb{Q,PH) ® OV^ph) 




mt 
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Compared to the fixed-order case in eq. (2.2), where we had ~ Aqqq, we now have 
p^ ~ which not only includes 6-quarks but also collinear gluons of that virtuality. 
Thus, as indicated in eq. (2.17), it is always an expansion in . In particular, matrix 

elements with external 6 quarks (e.g. in the matching calculation below) are expanded in 
the 772-6 0 limit. The coefficients Ci^b{Q, pn) contain the full dependence on the physical 

scale Q but are independent of the low scale mb- On the other hand, the operators do 
contain an implicit rub dependence since they involve massive 6-quark helds. 

In principle, one is free to reabsorb some of the neglected 0{mflQ‘^) corrections in 
eq. (2.17) into the coefficients. This corresponds to including some subleading power 
(subleading twist) corrections in the leading-power (leading-twist) result and letting the 
leading-power resummation act on them. However, we stress that to correctly include the 
subleading power corrections in the resummation requires extending the factorization in 
eq. (2.17) to the subleading order. As mentioned before, the power corrections in 
can be important in practice and should be added back such that in the fixed-order limit 
Pm —^ Ph we recover the fixed-order result of the previous subsection. This is discussed in 
detail in section 2.4. In the rest of this section, we do not indicate the power corrections 
for ease of notation. 

After the matching at pn in eq. (2.17), we want to evolve the theory from pn down to 
Pm- The renormalization of the PDF operators again gives rise to their RGE, the solution 
of which is given by DGLAP evolution, relating the operators at different scales, 

of ip) = Ulf {p, po) ® of {po) + Uf {p, po) ® of (/2o), 

0f{p) = uf{p,po) ® Of{po) + uf{p,po) ® 0f{p ,). ( 2 . 18 ) 

f5l 

The difference to eq. (2.5) is that now Uf = 5 and O^ contributes to the evolution, which 
we have written out explicitly. Taking the proton matrix elements on both sides yields the 
corresponding evolution of the PDFs from pQ to p in the theory above pm- Equivalently, 
we can use eq. (2.18) to evolve the Wilson coefficients from pn down to pm, 

Cj{Q,Pm) = Ci{Q,pH) ® Ulf{pH,Pm) + Cb{Q,PH) ® Ulf{pH,Pm) , 

CbiQ,Pm) = Ci{Q,pH) ® Uf{pH,Pm) + CbiQ,PH) ® Ull\pH,Pm) ■ (2-19) 

Next, at the scale pm, the operators of{pm) and of{pm) are matched onto the set of 
operators of{pm), which are precisely the ones appearing in eq. (2.2) and do not include 
a 6-quark operator, 

of{pm) = .Mjkimb,Pm) ® of{pm), (2.20) 

ofiPm) = -M-bkimb, Pm) ® of{pm) ■ (2.21) 

f5l 

By integrating out the 6 quark, the mb dependence implicit in the Oj b{pm) is now fully 
contained in the matching coefficients .Mjk{mb, Pm) and -Mbkimb, Pm)- In particular, the 
Ob operator does not exist in the theory below pm and its effects are moved into the Mbj 
coefficient. In addition, secondary 6-quark loops are integrated out, which corresponds to 
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switching the UV renormalization scheme for b quarks at the scale //m from MS to the 
decoupling scheme, and the Adij and Aibj contain the associated matching (threshold) 
corrections. 

The remaining steps now proceed as in the previous subsection. The operators (or 
PDFs) at fim still contain logarithms ln(^m/AQCD), which are resummed by using eqs. (2.5) 
and (2.6) to evolve them from up to fj,m- Equivalently, we can think of evolving the 
products Cx{Q, Hm)-M.xk{rnb, IJ^m) from /im further down to (with x = j,b). The final 
expression for the DIS cross section is then given by 


d^resum(g,^^) 

+ ( 8 ) +CbiQ,fiH) ® 

® ^kl ^a) ® /f' (pa) ■ 


® l^m) 

0 Mbk{mb, Hm)'^ 

( 2 . 22 ) 


The full cross section on the left-hand side contains large logarithms \n{Q/mb) and ln(mft/AQCD)- 
On the right-hand side these are factorized into logarithms ln{Q/ij.H) and liii{mb/pm), 
which are considered small and reside in the coefficients Ci^biQ-, Ph) and M.{mb, Pm), large 
logarithms Pm) and Pa), which are resummed into the evolution factors 

U^^\ph, Pm) and Pa), and finally logarithms ln(/rA/AQCD)5 which are absorbed 

into the PDFs at //a- We will refer to eq. (2.22) as the resummed (“resum”) result, since 
it has all logarithms ln(Q/mb) resummed. 

In the traditional 5F scheme, the resummed result in eq. (2.22) is written as 

da^^{Q,mb) = Cb{Q,PH) ® ff\pH,mb) + Ci{Q,pH) ® ff\pH,mb ), (2.23) 


where the combinations 


ft\rnb,PH) 

ff\mb,pH) 


^bj ^PH, pm) ® jkip^b, pm) {PH, pm) ® Pm) 

^ij ^PH, pm) ® ■f^jkijblb, pm) ~\~ U{^PH, pm) ® ■M-bkijblb, pm) 


® fk^ (Pm) , 

® fk ^ {Pm) , 

(2.24) 


f5l 

are interpreted as the evolved 5F PDFs including a PDF for the bottom quark /^ . To all 
orders in ag, eqs. (2.23) and (2.24) are simply a different way to write eq. (2.22). In practice, 
however, the evolution and matching corrections are always carried out to a certain finite 
order, where the different interpretations lead to different perturbative countings yielding 
different results. This is discussed in detail in section 2.3. Basically, in eq. (2.23) the 
5F PDFs are traditionallly regarded as external 0(1) inputs, and the perturbative order 
counting in ag is only applied to the coefficients Ci and Cb- In contrast, in eq. (2.22), we 
only regard the as external 0(1) quantities, while the perturbative order counting 

is applied to all terms in curly brackets. As we will see, one advantage of doing so is that 
this renders the order counting consistent between the resummed and fixed-order results, 
which facilitates their combination, as discussed in detail in section 2.4. 
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Figure 4. Schematic leading-order matching at fiu for heavy-quark production in DIS for mb ^ Q. 



Figure 5. Schematic NLO matching at for heavy-quark production in DIS for mb Q. 


We also note that the publicly available 5F PDF sets are constructed as in eq. (2.24), 
with the notable difference that the matching scale is commonly identified with and 
fixed to the heavy-quark mass, fj,m = rrih- However, it is clear from our discussion that 
IXm is a (in principle arbitrary) perturbative matching scale and it is important to keep it 
conceptually distinct from the parametric mb dependence. In our results, we will utilize 
the dependence to estimate the intrinsic resummation uncertainties. 

2.2.1 Matching at hh ~ Q 

The Wilson coefficients and Cb{Q,HH) are computed in perturbation theory by 

taking partonic matrix elements of both sides of eq. (2.17), 

{b\OmsiQ,mb)\b) = Cb{Q,fiH) ® {b\of\^iH)\b) + ® {b\of\^iH)\b), (2.25) 

{j\ODis{Q,mb)\j) = Cb{Q,HH) ® {j\of\nH)\j) + Ci{Q,HH) ® {j\of\^lH)\j) ■ (2.26) 

The calculation proceeds analogous to section 2.1.1. The essential difference is that now 
h quarks are present in the theory below and so we also have to consider external 
6-quark states to determine the Cb matching coefficient. As discussed earlier, the full- 
theory matrix elements on the left-hand side are expanded to leading order in vn^/C^. The 
partonic matrix elements on the right-hand side now lead to partonic PDFs similar to those 
of eq. (2.10), now including also 6-quarks, 

= {b\of\nH)\b) , flfbi'rnb^m) = {b\of\nH)\b) , 

flfjimb, fJ^H) = {j\of\fiH)\j) , flJ.{mb,HH) = {j\0^b\f^H)\j) ■ (2.27) 

To perform the matching, we expand both sides of eqs. (2.25) and (2.26) in powers 
of where all the pieces are expanded analogously to eq. (2.11). The 

leading-order matching is illustrated in figure 4. At LO, the partonic bottom PDF is 

fb^b^\z,mb,HH) = S{1-z), (2.28) 
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so the LO bottom-quark coefficient is directly given by the LO diagram on the left of 
figure 4, 

= lim {b\Oms{Q,mb)\bf^ . (2.29) 

The limit mb —)• 0 explicitly highlights that the full-theory matrix element is expanded in 

The NLO matching is illustrated schematically in figure 5. At this order, there are 
1-loop and real-emission corrections to the bottom-quark LO contribution as well as a 
contribution from a gluon channel. The partonic PDFs at NLO for finite mb are (see 
e.g. ref. [52]) 

fl/b^\^,mb,fiH) = 2 CfO{z) 

2 

fb/g^\^^^b, m) =2TFe{z)Pqg{z)\\l^, (2.30) 

with 


1 + z^ 
1-z 


In - 




m?(l — zY 


- 1 


J + 




(2.31) 


Together with the LO 6-quark coefficient in eq. (2.29), the NLO matching coefficients are 




(2.32) 


The logarithms of mb inside the matrix elements of the DIS operator precisely match those 
in the partonic PDFs in eq. (2.30), such that the mb —>• 0 limit is finite. The reason is 
that for the matching at ///f, the finite bottom mass is nothing but an IR regulator for the 
collinear divergences associated with bottom quarks, which cancels in the matching. 

Since the matching coefficients are independent of the IR regulator, we can also take the 
mb —?■ 0 limit at the beginning, as long as we use another IR regulator, such as dimensional 
regularization. In this case, the computation of the coefficients Ci^b becomes much simpler 
since there is one less scale involved. The partonic PDFs are then the usual ones in pure 
dimensional regularization, completely analogous to eq. (2.14), 


ft}i^\z) = -\2CFe{z)Pqq{z), 
ft}l^\^) = -\‘iTFe{z)Pgg{z), 


with 


Pqq{z) = 


Oil - Z) 


1 + Z^ 

1 - z 


and Pqg{z) as in eq. (2.31). The NLO coefficients are then given by 


= {b\OBis{Q.^)W -criQ^b^H)®fb,b 

and are precisely the same as in eq. (2.32). 




do). 


[5](1) 


(2.33) 

(2.34) 


(2.35) 
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Figure 6. NLO matching at 


2.2.2 Matching at ~ 

To compute the matching coefficients Mij at the low scale we calculate matrix elements 
of both sides of eqs. (2.20) and (2.21) with the same external partonic states. Using the 
definitions in eq. (2.10) and eq. (2.27), we have 

/|/|j(^6)/^m) — IJ-rn) ® /jy^fc(l^m) ? 

IJ-m) = Mbj{mb, Urn) ® /]/fc(/^m) ■ (2.36) 


Now, the b quark cannot appear anymore as an external state, since it is integrated out on 

rci rci 

the right-hand side. (Hence, there are no equivalent matching equations for or f)^l-) 
At the same time, nib is now the hard scale which appears in the matching coefficients 
(and cannot be set to zero). The matching coefficients My are known fully to 0{ag) [53]. 
(They are also known partially to O(a^), see e.g. refs. [54, 55] and references therein.) 
Expanding eq. (2.36), the LO matching is simply 

M£)(z) = M(J(z) = 5(1 - z), Mg) = Mg) = Mg) = Mg) = 0. (2.37) 


At NLO, there are nontrivial matching conditions for Mg) and Mg), which are illustrated 
in figure 6, 

Mg)(2;, nib, Hm) = mb, Hm) - 4/gg^) = In ^ '^(1 “ > 

2 

Ml^g\z,mb,fXm) = fljg^\z,mb,fim) = 2TFe{z) Pqg{z) In . (2.38) 

Note that the precise number of flavors used in ag here is an O(a^) effect, which will then 

(2) (2) 

also generate nontrivial matching conditions for Mgg and -^bg ■ 


2.3 Perturbative expansion and order counting 

We now discuss the perturbative counting for the cross section for the two scale hierarchies 
in figure 1. Since the gluon and light quark PDFs at the scale are nonperturbative 
objects fitted from data, we make the standard assumption and count them as external 
0(1) quantities, 

/gkl^A) ~/gk/^A) ~/g)(/^A) ~ 0(1), q = d,u,s,c. (2.39) 
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To determine the cross section to a certain perturbative accuracy, a perturbative count¬ 
ing should be applied to all remaining terms in the cross section that are computed in 
perturbation theory. In the fixed-order case ~ Q, this implies that the standard per¬ 
turbative counting in terms of powers of ag appearing in the hard matching coefficients 
applies. On the other hand, for mt <C Q, the perturbative counting should be applied also 
to the matching coefficients at and the evolution factors between and Hm- We will 
argue that for phenomenologically relevant hard scales this implies that an appropriate 
perturbative counting takes the effective bottom PDF to be an 0{as) object. 

2.3.1 Fixed order 

First, recall that the DGLAP evolution factors ^ 2 ) resum single logarithms of the 

ratio li \/[12 to all orders in a^. For this purpose, one performs a logarithmic counting 
where one expands in powers of while counting as\'a{^i/^ 2 ) = oigL ~ 1. That is, one 
formally counts L ~ 1 /ug ■ We can then write 

UifMi,fi2) = U^^iagL) + a, U^'^^{agL) + U^^'^^{agL) + ■■■ 

~ 0(1) + 0{ag) + 0{al) +•••, (2.40) 

where are functions of agL to all orders in ag. Combining this with eq. (2.39), we 

can also count the evolved 4F PDFs as 0(1) quantities 

® ~ 0(1). (2.41) 

The PDFs evolved at N^LL are then usually called N^LO PDFs. Note that while the 
evolution mixes the PDFs, it does not induce a parametric difference between the light- 
parton PDFs. Also, in the limit —)■ ^2 we have 6ij. Hence, we can generically 

treat as external 0(1) quantities for any /r regardless of how large the logarithms 

IJ. 2 ) actually are, and this is the standard praxis. 

For the fixed-order (4F) cross section in eq. (2.8), the perturbative counting in ag 
is then directly applied to the Di coefficients, so it has the perturbative expansion [with 
ag = (ijh )/as in eq. (2.11)] 

LO (FO, 4F) da^^{Q,mb)= ag ® 

NLO (FO, 4F) + of ^(Q, m^, fin) ® 

NNLO (FO, 4F) + 4 Df'\Q,mb, fin) 

+ •••. (2.42) 

Taking into account eq. (2.40), obtaining the cross section at an accuracy of order 
(N^LO) then requires the N^LO matching coefficient with the N^LL evolution (i.e. N^LO 
PDFs). Here the order is counted relative to the lowest nonvanishing order, which for 
heavy-quark production in DIS is 0{ag). 

When combining the expansion in eq. (2.40) with the a* expansion of Di{Q,mb, hh), 
one could in principle reexpand the product of the two series. In practice, this is usually 
not done, since the /]^^(/^h) are treated as external 0(1) inputs as mentioned above. 
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2.3.2 Resummation 

We now discuss the perturbative counting for the resummed cross section in eq. (2.22). 
First, we can use the same arguments as in eq. (2.41) to treat the 4F PDFs at as 0(1) 
inputs. We then have to consider the perturbative counting for the terms in curly brackets 
in eq. (2.22). The situation is more subtle now due to the presence of the additional scale 
IXm ~ ^b- Depending on the hierarchy between hh and Hm-, there are two different options 
of how to count the evolution factors , ^J-m)■ 

• For very large hierarchies, we can use a strict logarithmic counting, in which case 
otgL 1 and eq. (2.40) generically applies to all the evolution kernels so 

(2-43) 


• For intermediate hierarchies /im ^ fJ-H, the resummation can still be important, so we 

still use eq. (2.40) to organize the logarithmic order of the resummation. However, 

we should also take into account that in the limit fJ-H the off-diagonal mixing 

evolution kernels vanish 0 and similarly for This is 

because their fixed-order expansion starts at order agL rather than 1, so they are 

[5] [5] 

suppressed by an overall factor of agL relative to the diagonal and Ugg'. Therefore 
we count 

Utj ~ ~ • (2.44) 


The counting in eq. (2.43) corresponds to the traditional 5F scheme. With this counting 
and using eqs. (2.37) and (2.38), the evolved 5F PDFs in eq. (2.24) have the perturbative 
expansion 


#(7716,/in) = 


l-H 


) /^m) T 


47r 


^gb ) l^m) ® '^bg T 




0(1) + Oiag), 

0 ( 1 ) 


47r 


1 l^m) ® "^bg i^b^ l^m) T 




0{ag). 


} ® #(Mm) 

} ® #(Mm) 

(2.45) 


Hence, they are treated as external 0(1) quantities. The resummed result is then written 

f5l 

as in eq. (2.23) and has the perturbative expansion [with an = (/i77)/(47r)] 


LO (5F) d(T{Q,mf,) 
NLO (5F) 

NNLO (5F) 


® fnrnb,f^H) 


[5]/ 


+ an 


(8) #(r?l6,^77) -h ® #(7716,/ill) 


45] 


(1)7 


45]/ 


+ (8) /f'(m6,/ij7) -4 cf\Q,iXH) ® 


+ 


(2.46) 
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The N^LO cross section then requires using the N^LO 5F PDFs, which are given by the 
expansion of eq. (2.45) to 0{ag) together with the N^LL evolution factors. 

A rough numerical estimate shows that for ~ nT-b ~ 5 GeV the first case eq. (2.43) 
applies for hard scales fj-n ^ 1 TeV. Thus, the second case in eq. (2.44) is more appropriate 
for our purposes. This is also confirmed by the fact that for fin ~ 0(100 GeV) and standard 
PDF sets one finds that numerically /]^^(^h) ^ Adopting this counting, the 

resummed result in eq. (2.22) has the perturbative expansion 


LL (resum) 
NLL (resum) 



Here, an = ai^^(^_ff)/(47r) and am = cJs\fJ-m)/, and for notational simplicity we have 
suppressed the convolution symbols and all arguments (which are as in eq. (2.22)). In the 
contributions proportional to we have also counted Ugq ~ ag and Ubq ~ a^. Note 
that, in the region where this counting applies, an and am can be regarded as being 
parametrically (and practically) of the same size. 

From eq. (2.47) we see that the counting in eq. (2.44) leads us to include the matching 
terms Cg^'^ and , which provide the boundary conditions for the RGE, already at 
the lowest order, i.e. one order lower compared to the 5F. Furthermore, any cross terms 
in eq. (2.47) from the matching at hh and Hm are expanded against each other. In other 
words, compared to eq. (2.46), we do not have overall 5F PDFs, but rather the contributions 

f5l 

~ Ulj IJ-m) ® -Mjkiiam) making up the 5F PDFs in eq. (2.45) are expanded together 
with the hard matching coefficients. As we will see in the next subsection, these features 
enable us to have an easy and smooth transition to the fixed-order result. Note that this 
is quite similar to how the primed resummation orders N^LD are implemented in the 
resummation for differential spectra, see e.g. refs. [56-59], where this facilitates a clean and 
smooth transition to the fixed-order result. 

We can of course collect the terms proportional to Cb and Ci in eq. (2.47) into effective 
PDFs, which we denoted as fb and fi to distinguish them from the standard 5F PDFs in 
eq. (2.45). With the counting in eq. (2.44) we then have 


Mmb,HH) = (8) /f'(^m) + ' ' ' 

0(1) + 0{as), (2.48) 

fb{mb,flH) = ulf {fJ-H, {rUb, Urn) ® 4‘^k/^m) H- 

0{as) + 0(as) +0{al). 
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Thus, in the region of scales we consider, the effective 6-quark PDF should be treated as 
an 0(as) object, while the gluon PDF still starts at 0(1). Note though that the 0(as) 
terms in fg are different from those in fg^^. For example, the <8) term in eq. (2.45) 

counts as 0(as) in fg^^ while it only appears at 0(ag) in fg. Since the light-to-light Mgg 

and the light-to-heavy Mhg matching functions are needed at different relative orders in 

fbi'mb, fJ'H), this definition of the effective bottom PDF differs with respect to the usual 
f5l 

fb \^^H) in eq. (2.45). Hence, in our numerical implementation we cannot use the 6-quark 
PDF from the standard 5F PDF sets. Instead, we need to construct fbiiT^b, ourselves. 
We do so by creating PDF grids that have the matching coefficients at the required order but 
the same order in the evolution factors. The technical details are discussed in Appendix B. 

Denoting with the truncation of the effective PDF fi^b to 0{ag), we can write the 
NLL result in eq. (2.47) using eq. (2.48) in a compact form as 

da^^^{Q,mb)= anC^^HQ, i^ih) ^ fyHrrib, ^ih) + ^ m) 

+ ® fl^^{mb-, iJin) + aHC\^\Q-,^iH) ® fl^^{mb-, ij-h) - 

(2.49) 

Here, we still consistently drop any higher-order 0{ag) cross terms in the product of 
coefficients and effective PDFs by keeping the effective PDFs to different orders in the 
different terms. As already mentioned, this is important to ensure a smooth transition to 
the fixed-order result in the limit Hm ~ fJ-H- 

On the other hand, for the bbH hadron collider process we are eventually interested in 
in section 4, we will have two PDFs and the practical implementation of the strict expansion 
gets quite involved. Therefore, as long as we are only interested in the phenomenologically 
relevant region ~ mb hh ~ mn, we can also keep the higher-order cross terms to 
simplify the practical implementation. We then have 

da~{Q,mb) 

= 'Si fg{mb,^J,H) + fJ-n) S fb{mb, ^J-H) 

+ aHCl‘^\Q, fJ-H) S MmbjfJ-n) + aHCl^\Q, fJ-n) S fb{mb,fJ-H) 

+ •••. (2.50) 

Once we allow keeping higher-order terms, we can also further simplify the practical im- 

[51 

plementation by replacing the effective PDFs fi^b above by standard 5F PDFs fi^. These 
must then be of sufficiently high order such that they include all necessary matching cor¬ 
rections as required by our perturbative counting. However, we note that whenever one 
keeps higher-order terms for practical convenience, one should check that this does not 
have a large numerical influence on the results in the kinematic region of interest. We will 
come back to this in section 4. 

We stress, that even when keeping higher-order cross terms, the perturbative counting 
is still performed for both Ci^b and fi^b with fb counted as 0{as)- So even though the leading 
term in (7^°^ is 0{a^), the resummed result starts at 0{as)- Comparing to eq. (2.42), the 
resummed result has a perturbative counting consistent with the fixed-order result. It 


LL (resum) 
NLL (resum) 
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precisely corresponds to a resummed version of the fixed-order result in the mb —>• 0 limit. 
This organization and implementation of the resummation is one of the main ways in which 
our approach differs with other approaches. This will be discussed further in section 3. 


2.4 Combination of resummation and fixed order 

In sections 2.1 and 2.2 we have derived results for the heavy-quark production cross section 
in DIS that are relevant for two different parametric scale hierarchies. The fixed-order result 
in eq. (2.8) is relevant for mb ~ Q, as it keeps the exact mb dependence to a given 
fixed order in Og, but does not include the all-order resummation of logarithms ln{mb/Q). 
The resummed result in eq. (2.22) is relevant for mb Q, as it resums the logarithms 
ln{mb/Q) to all orders in ag, but neglects any mb IQ power corrections that vanish for 
mb —7- 0. These two results represent two ways of computing the same cross section. In this 
section, we combine these two results and obtain our final result accurate for any value of 
rriblQ. 

We follow the usual approach for combining a higher-order resummation with its cor¬ 
responding fixed-order result. We write the full result for the cross section as 

da = da’’"™“ + da“°“". (2.51) 


Here, the nonsingular cross section dcr'^°“® contains all contributions that are suppressed 
by 0{mbjQ) relative to da’’®™™ and vanishes in the limit mb —)• 0. With this condition, da 
automatically contains the correct resummation in the mb —)• 0 limit. 

Furthermore, we require that the fixed-order expansion of eq. (2.51) reproduces the 
correct FO result, including the full mb dependence. Therefore, 


da"°^'^ = da^° - da“°s , 




daresum 


FO ’ 


(2.52) 


where the singular contributions da®™® are obtained from the fixed-order expansion of the 
resummed result to the desired order in Og. For da“°°® to indeed be nonsingular and 
vanish for mb —)• 0, da®™® must contain all singular contributions in da^®, i.e. all terms 
that do not vanish as mb —)• 0. This in turn requires that the resummation to a given order 
fully incorporates all these fixed-order singular terms. In this sense, the resummed result 
should be consistent with the fixed-order result. This condition is precisely satisfied by our 
resummed result with the perturbative counting used in eqs. (2.47) and (2.50), for which 
the (N)LL result contains the full (N)LO singular terms, as we will see below. 

To explicitly identify the nonsingular terms, we need a meaningful and consistent 
comparison between da^® and da®™®, which means we have to write both in terms of the 
same external 4F PDFs and expand both in terms of the same Og. For this purpose, it 
is most convenient to use /]^^(/.t_H') as in eq. (2.42) but perform the expansion in terms of 
as in the resummed result eqs. (2.47) and (2.50). First, for da^®, we can simply 
change the 6-quark renormalization scheme for Ug used to computed the Di matching 
coefficients in eq. (2.42) from the decoupling scheme to the MS scheme. This leads to 
modified Wilson coefficients Di{Q, mb, fJ-n) —^ mb, ^ih) which are now expanded in 

terms of the same «« (//_«■) as is used in Ci{Q, ^H),Cb{Q, ^h)- From the point of view of 
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the fixed-order calculation, this is actually the more appropriate expansion for mb < hh- 
Next, the singular cross section can be easily obtained by evaluating the resummed result 
in eq. (2.47) or eq. (2.49) at 


da'^^s = da’’ 


IFO 


= da’’ 




Cj{Q,^iH)'SiMji{mb,^iH) + Cb{Q,lJ.H)'SiMbi{mb,iiH) (2.53) 


Finally, the fixed-order nonsingular cross section is given by 
= dcT^° - da"“s 

= Df^{Q,mb,iJ.H) - Cj{Q, hh) <Si Mji{mb, hh) - Cb{Q, fJ^n) 'S' Mbi{rnb, hh) 

(2.54) 


At each order in a^, all singular terms in LiNS exactly cancelled by the corresponding 
singular terms from the resummed result, such that d(T°°'^® is free of collinear logarithms 
and vanishes as mb —>■ 0. 

We stress that the statement d(T’’®®““|FO = utilized above is quite non¬ 

trivial and crucially relies on the fact that with our perturbative counting in the resummed 
result all the matching corrections Mij are always included to sufficiently high order (ba¬ 
sically to the same order in as to which we have to expand the evolution kernels) such 
that the fim dependence precisely cancels in da®™® to the given order in Og to which we 
expand. Once we know that this is the case, we can pick any we like to perform the 
fixed-order expansion of da^®®“™. The choice Hm = is then the most convenient, since 
all the evolution kernels become trivial. For example, at LL we have 


^\g {fJ-H J (nib, fJ^rn) 


= anM^f^^mb, ^jlh) 

fJ'Tn — 


(2.55) 


and therefore 


(i^singLO ^ 


LL| 


LO 




= CLH 


c^^HQ, m) + c^°\q, m) s M[l\mb, i^h)] s 4‘^k/^H), 


da 


nons LO 


= an 


(Q, mb, ^ih) - {Q, hh) - O'® {Q, hh) S {mb, 


iW, 


bg 




(2.56) 


Comparing to the matching conditions in eqs. (2.32) and (2.38), we can see explicitly that 
the last two terms in square brackets in da“°°® precisely reproduce the singular mb —)• 0 
contributions of D^^^^\Q,mb, Similarly, at NLL we have 

daS™gNLO = daNLL| . (2.57) 

Note that in eqs. (2.56) and (2.57) we have implicitly assumed that the resummed result is 
taken as in eqs. (2.47) and (2.49), with all cross terms consistently expanded. Otherwise, 
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e.g. when using eq. (2.50), any higher-order cross terms then need to be dropped at the 
level of d(T®™® to avoid introducing spurious uncancelled singular terms in dcr°°'^®. 

So far, the nonsingular corrections are expressed in terms of 4F PDFs at the hard scale 
Hh, while the resummed cross section is necessarily written in terms of 4F PDFs at fim or 
effective as in eq. (2.49) or eq. (2.50). To simplify the practical implementation 

it is desirable to only deal with a single set of PDFs. For this purpose, we can choose to 
write the nonsingular contributions in terms of only light-parton effective PDFs /i(mf,, ^h) 
as 

d^nons ^ AC'f°“"(Q, mb, fin) ® fi{mb, i^h) , (2.58) 

where the new coefficients ACf°'^^{Q,mb, fj-n) are fixed by equating this to eq. (2.54) at 
each order in Og- This has a unique solution, since the nonsingular contributions are by 
definition a FO contribution, so at the terms in eq. (2.54) always have the form 

_ (jM _ ... j (g) Therefore, we can naturally associate them with the gluon and 
light-quark PDFs /j. We can then absorb the nonsingular corrections into the light-parton 
coefficient functions by taking 

CiiQ, ^ih) Ci{Q, mb, ^ih) = Ci{Q, fin) + AC'f°“‘'(Q, mb, fJ-n) , (2.59) 

while keeping the Cb coefficient unchanged. Equivalently, we can replace C, —?• Ci every¬ 
where and impose the condition 

Ci{Q,mb, yLu) ® Mij{mb, ijlh) = Df'^{Q,mb,^lH) - Cb{Q, iih) ® Mbj{mb, ^ih) , (2.60) 

such that eq. (2.54) vanishes. 

The above shows that we can choose to absorb the nonsingular contributions into the 
resummed result by modifiying the matching coefficients at fin- The condition in eq. (2.60) 
implies that the light-parton coefficients Ci{Q, mb, fJ-n) can be obtained from the matching 
at fj,H in section 2.2.1 without taking the mb —?■ 0 limit in the light-parton full-theory 
matrix elements, while for all bottom contributions and coefficients the mb —?• 0 limit is 
still taken. 

We can now write the final result for the cross section as 
du = da’'"™” + da““" 

= Ci{Q, mb, iih) ® fi{mb, hh) + CbiQ, Hh) ® fb{mb, fJ-n) , (2.61) 

which now uses the effective PDFs fi^b throughout whilst capturing the full nonsingular 
corrections. The same perturbative counting as in eqs. (2.49) and (2.50) still applies, which 
now gives 


da{Q,mb) 

LO-hLL = aHC^^'>{Q,mb,^J.H) ® fg{mb,fiH) + ^J-H) ® fb{mb, fin) 

NLO-hNLL ^ {Q, mb, hh) ® fi{mb, hh) + (Q, I^h) ® fb{mb, fJ-n) 

+ •••, ( 2 . 62 ) 
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where the choice of expanding the cross terms or not is kept implicit and will determine to 
what order the PDFs are kept in each term. We emphasise that in this form the result is 
very convenient to implement, since it essentially only requires the fixed-order result (after 
changing the 6-quark renormalization scheme for and the massless resummed result. 
In section 4 we will apply this strategy to the bbH cross section and provide further details 
on the construction of the coefficient functions. 

By choosing to absorb dcr°°'^® into the matching coefficients in the resummed result, we 
effectively let the leading-power resummation also act on the nonsingular corrections. This 
introduces power-suppressed higher-order logarithmic terms, which however are beyond 
the order we are working at. In particular, this does not include the correct resummation 
of power-suppressed logarithmic terms. (This would require the extension of do-''®®'!™ to 
subleading order in rriblQ, which is well beyond the scope of this work, and also very 
likely irrelevant at the current precision.) Fundamentally, we only have control over the 
nonsingular corrections at the level of their hxed-order expansion. The above procedure 
to include the nonsingular contributions is not unique, and while physically motivated, is 
ultimately driven by practical convenience. We would like to underline that alternative 
choices are in principle possible, provided they do not change the resummation in the 
mb —)• 0 limit and reproduce the correct fixed-order expansion, in which case they will 
effectively differ by power-suppressed higher-order logarithmic terms. ^ We will come back 
to this point in section 3. 

From the discussion so far, it is clear that transition between and dcr^® is 

controlled by the scale fim- To provide a smooth transition between the resummation 
and fixed-order regions, this scale is promoted to a mb-dependent profile scale —>• 
IJ'mirnb, fJ'H)■ It has the properties that in the resummation region for mb <C Q it has 
the canonical resummation scaling ~ w-t, while in the fixed-order region mb ~ Q it 
approaches Hm —^ fJ-H, such that the resummation is turned off there and the fixed-order 
result is recovered, with a smooth transition in between. The fact that it is possible to 
control this transition between limits with a single scale, makes our predictions in the 
transition region robust and, moreover, variation of this scale and of its functional form 
provides a solid handle on the associated theoretical uncertainties. The precise definition 
and variations of the profile function are discussed in detail in section 4.3 for the case of 
bbH production. 

3 Comparison to existing approaches 

In the previous section we used a systematic field-theory analysis to derive a result for 
the heavy-quark production cross section in DIS accurate for all possible scale hierarchies 
from m <C <5 to m ~ Q. Various approaches to the same problem are available in the 
literature, which go under the name of variable flavor number schemes (VFNSs). In this 

^In general, one could write the nonsingular contribution in terms of both light-parton and bottom 
PDFs, fiHmb, and in this case there would not be a unique solution to eq. (2.58). This gives rise 
to several (equivalent) possibilities of writing the final result, and this generates some of the differences 
between the various VFNSs. 
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section, we briefly compare the existing schemes to our result. In what follows, we do not 
attempt to give an in-depth review of the different schemes but rather focus on similarities 
and differences with respect to our EFT result. For reviews of the different schemes in the 
literature we refer to refs. [60, 61] and sect. 22 of ref. [62]. 

VFNSs can in principle differ in various aspects. The first is the general construction, 
namely how resummation of collinear logarithms is achieved and how nonsingular power 
corrections are included. Secondly, they can differ in how the perturbative counting is 
performed, that is, which of the various perturbative ingredients are included at a given 
order. Third, they can differ in how the heavy-quark threshold is implemented, which 
in our language corresponds to the exact choice of the low matching scale [Xm- The first 
aspect is the one that primarily distinguishes the different schemes, while the remaining 
two aspects are more related to choices made within each scheme. Here, we compare to 
the choices often used in the literature. We stress though that these choices correspond 
to how a particular scheme has been used or implemented in practice, but (in most cases) 
they do not necessarily represent restrictions of a particular scheme itself. 

3.1 Construction 

We start by discussing the differences in the basic construction of the cross sections. For 
mi) > Q (“below threshold”) all schemes use the same fixed-order 4F result in eq. (2.8). For 
mi) < Q (“above threshold”) the various schemes construct their cross sections as follows: 

• Zero-Mass (ZM). In this approach, the massless resummed result in eq. (2.23) is 
used for m^ < Q, while nonsingular power corrections are neglected at any order in 
as- Hence, this scheme is only expected to be accurate for mi, <C Q. Since power- 
suppressed contributions are not included, it is not accurate close to the heavy-quark 
threshold and does not reproduce the full fixed-order result. For this reason, we do 
not discuss it further. 

• ACOT [1, 2, 5]. The ACOT scheme is based on the idea that the power corrections 
can be fully included in DIS at the level of the matching at the hard scale 

eq. (2.17), by generalizing it such that power corrections in mhjQ are included in 
the definition of the Wilson coefficients. The heavy quark is considered as an active 
flavor and the quark mass dependence is retained at each matching step, yielding 

do- = Ci{Q,mb,^XH) ® fl^\mi„fXH) + Cb{Q,mi), hh) ® fl^\rni), hh) • (3.1) 

The Ci^i){Q, mi), ^h) incorporate the nonsingular contributions and reduce to the orig¬ 
inal Ci^biQ^ (J-h) in the ruft —)• 0 limit. In contrast to eq. (2.61), the heavy-quark con¬ 
tributions in eq. (3.1) are computed with a massive on-shell heavy quark in the initial 
state. To account for the massive kinematics including the presence of a massive (un¬ 
resolved) heavy quark in the final state, the heavy-quark Bjorken-x can be rescaled, 
leading to a variant of this scheme called ACOT-x [7, 63]. While the validity of 
ACOT in DIS can be based on including heavy-quark masses in the hard-scattering 
factorization [5] , its extension to the case of two incoming hadrons is problematic due 
to the massive kinematics, see Appendix C. 
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• S-ACOT [6]. The fact that is not independent of (see eq. (2.24)) allows one 

to move power corrections between Cf, and Ci without spoiling the formal accuracy of 
eq. (3.1) [5, 6]. This was used to construct a simplified variant of ACOT, in which the 
heavy-quark Wilson coefficients are computed in the massless limit, Ch{Q, mi,, —)• 

CbiQ, fJ-n), while the full mass dependence is retained in the (modified) light-parton 
coefficients. This is evidently equivalent to how we include the nonsingular corrections 
in eqs. (2.60) and (2.61) for practical purposes. To account for the massive heavy- 
quark kinematics, a x-rescaling is also applied, leading to S-ACOT-x [7, 13].® In 
ref. [31], a modification of ACOT, dubbed m-ACOT, is used for the case of two 
incoming hadrons, where the massless limit is applied only to channels with two 
incoming heavy quarks, while the mass dependence is kept in heavy-light and light- 
light channels. 

• TR [3, 8]. The TR scheme is defined by requiring that the fixed-order result, after 
being expressed in terms of 5F PDFs, corresponds to the resummed result up to 
power-suppressed contributions. This requirement fixes the singular contributions. 
However, there is still freedom for the treatment of nonsingular terms, and this is 
fixed by making a choice such that the coefficient functions obey a sensible threshold 
limit. The result is hence different from both ACOT and S-ACOT. Due to the choice 
of perturbative counting, a discontinuity exists at threshold which is removed by 
adding a Q-independent contribution to the result above threshold. Though this 
contribution is formally higher order, it can be sizeable, even far from threshold. 
The presence of the constant terms complicates the generalization of this scheme to 
higher-orders and to hadron-hadron collisions. 

• FONLL [11, 32]. This scheme is constructed by adding the massless resummed result 
to the full fixed-order result and consistently subtracting the double counting order- 
by-order in as- The fixed-order contribution is rewritten in terms of 5F PDFs with 
the resulting ambiguity fixed through the choice that only light channels contribute, 
as we have also done in section 2.4. The double-counting terms are equivalent to the 
singular terms in our notation, and remove from the fixed-order result its massless 
limit, i.e. all its terms that do not vanish in the mi, —)■ 0 limit. The FONLL procedure 
is thus equivalent to adding the d(T“°“® to the resummed result. This also makes 
the FONLL construction formally equivalent to S-ACOT. Finally, a damping factor, 
which performs the same function as the x rescaling in S-ACOT, is used to suppress 
higher-order spurious contributions and guarantee continuity at threshold. 

From the point of view of the all-order resummation, all these schemes are equivalent, 
as they all include the same resummation. As discussed in section 2.4, the minimal and 

®The X rescaling is not uniformly used in the literature. In some cases, the rescaling only takes into 
account the resolved b quark, whose kinematics is massive in ACOT but massless in S-ACOT (this variant 
only applies to S-ACOT), while in other cases, the x rescaling takes also into account the unresolved b 
quark, whose massive kinematics is not taken into account even in ACOT (this variant applies both to 
ACOT and S-ACOT). 
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formally correct result above threshold is given by eq. (2.51) as da = -)- dcr°°'^®. 

This result is formally correct in the sense that it correctly resums the collinear massive 
logarithms and correctly includes the full mass dependence and kinematics at fixed order. 
It is minimal in the sense that the nonsingular corrections d(T“°°® are unambiguous and 
unique when written in terms of 4F PDFs as in eq. (2.54), and are strictly included at fixed 
order, while the resummation strictly only includes leading-power terms. 

As discussed in section 2.4, there is an ambiguity when one tries to partially or fully 
absorb the nonsingular contribution into the resummed result, which amounts to expressing 
them in terms of 5F PDFs. The primary perturbative ingredients of ACOT, S-ACOT, TR, 
and FONLL are the same and they only differ in the way by which they fix this ambiguity. 
This ambiguity corresponds to power-suppressed higher-order logarithmic terms. Hence, 
these schemes can be regarded as formally equivalent up to such terms, which are beyond 
the considered formal accuracy. 

3.2 Combination of the ingredients 

We now move to the second source of scheme differences, namely how the perturbative 
counting is performed. By construction, the coefficient functions of (S-)ACOT, TR, and 
FONLL differ from each other and to those in our EFT result by formally higher-order 
contributions. Therefore, the largest differences between the approaches arise from the 
perturbative counting. In all practical implementations we are aware of, the perturbative 
counting used by each scheme is as follows: 

• ACOT-like schemes, used in the CTEQ family of PDE fits, construct perturbative 
expansions in the usual way by counting explicit powers of as in the coefficient 
functions. As a result, for DIS at LO (a^) the result below threshold is zero, while 
above threshold it is nonzero due to the heavy-quark initiated contributions (C^^^). 
At NLO, the gluon-initiated contribution (Cg^^) starts to contribute, as do the as 
corrections of the heavy-quark contributions 

• The TR scheme, used in MSTW and HERAPDE fits, is somewhat different as it com¬ 
bines the orders such that the lowest nonvanishing order below and above threshold 
appear at the same time. This means that at LO the result below threshold is the 
0(as) gluon-initiated contribution, while above threshold it is the O(a^) heavy-quark 
initiated contribution. The additional Q-independent term added above threshold is 
formally of higher order and does not affect the counting. 

• EONLL, used in NNPDE fits, also adopts the standard perturbative counting. The 
NLO and NNLO results are called EONLL-A and EONLL-C respectively. There is 
an intermediate result, EONLL-B, where the fixed-order terms are computed to order 
Qfg (NNLO) but the massless contribution is only included at order ag (NLO). 

None of the schemes discussed above adopts a perturbative counting which is directly com¬ 
parable to our approach of performing the counting on the full perturbative part of the 
cross section including both evolution and matching. In particular, it implies that the 
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Figure 7. Comparison of the construction of LO and NLO results in different counting schemes. 
Only representative diagrams at a given order and for a given channel are shown. Notice that, 
while the diagrams appearing in the Q < fim boxes contain collinear logarithms due to the heavy 
quark, the latter are subtracted in the diagrams appearing in the Q > boxes. Co,i are the 
Q-independent terms present in the TR-scheme that ensure continuity at threshold. We also point 
the reader to a very similar table in ref. [61]. 


effective heavy-quark PDF should be counted as an 0(as) object. Since the perturba¬ 
tive counting used by the different VFNSs summarized here does not distinguish between 
heavy-quark and light-parton initiated contributions, this difference in order counting is a 
principal difference in our approach. In figure 7 we summarize the perturbative counting 
adopted in the (S-)ACOT, TR and FONLL schemes as well as the counting we propose. 

As argued in section 2.3, our order counting is well justified theoretically and appro¬ 
priate for a wide range of scales (including scales appropriate for DIS experiments in the 
case of both bottom and charm quarks). It also has several advantages. As we will see 
in section 4, one of these is that the perturbative convergence tends to be improved with 
reduced uncertainties from the hard-scale matching. Another advantage, as highlighted in 
section 2.4, is that it facilitates a smooth transition to the fixed-order result at the heavy- 
quark threshold (provided the counting is strictly applied and higher-order cross terms are 
neglected), without the need of any rescaling or damping factors. 

3.3 Matching scale dependence 

Finally, the last important difference between the existing schemes and our approach is the 
position and treatment of the heavy-quark threshold. In all applications we are aware of. 
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the threshold is always set equal to the heavy quark mass, i.e. effectively the resummation 
scale Hm is fixed to Hm = 'mb- This is both the scale at which the low-scale matching 
is performed and also the scale at which one switches from the fixed-order result to the 
resummed result. 

Recently [61, 64], it has been suggested to consider an additional switching scale fis > 
fj,m, at which the computation switches from a fixed-order result to a resummed one, but 
nevertheless keeping the matching at a different (lower) scale. The effect of this choice is to 
delay the use of the resummed result, perhaps to a region where mass effects are negligible, 
though the transition between the resummation and fixed-order regions is not guaranteed 
to be smooth. 

In this work, we exploit the dependence on the matching scale to explicitly control 
the transition to the fixed-order result as rrn, —)• Q, and furthermore to estimate the intrinsic 
perturbative uncertainty in the resummation and matching procedure. This is in fact the 
standard practice in resummed calculations involving different resummation scales. This 
uncertainty should be taken into account as part of the total perturbative uncertainty in 
the result, which is typically not the case in existing approaches. 

4 Higgs production in association with b quarks 

In this section, we extend the framework presented in section 2 to hadron-hadron collisions 
and apply it to the bbH process, i.e. Higgs-boson production in association with b quarks. 
Specifically, this process can be defined as Higgs production via the bottom Yukawa cou¬ 
pling Yf,, with all other Yukawa couplings set to zero. (As discussed in section 4.2, we do 
not include the 6-quark loop contributions that are usually included in the gluon-fusion 
process. There we also comment on the inclusion of YtYb interference terms that are usually 
regarded as part of the bbH process.) 

The bbH process makes up only a tiny fraction, ~ 1%, of the total Higgs production 
cross section in the Standard Model (SM). It is nevertheless an interesting process within 
the SM, since the total bbH cross section is comparable to the total pp —)• ttH cross 
section for LHC energies, and because it provides direct access to the bottom Yukawa 
coupling. Furthermore, this process may be sensitive to new physics effects, since in many 
BSM scenarios, such as two-Higgs-doublet models with large tan/3, the Higgs coupling to 
bottom quarks can be enhanced. 

In the SM, the cross section in the massless 5FS is known at NLO [26, 27] and 
NNLO [28, 29], and in the 4FS at NLO [24, 25]. NLO predictions matched to parton show¬ 
ers for bbH production have been studied in both 4F and 5F schemes [65]. The 4FS and 
5FS calculations can lead to very different results, with cross sections differing by as much 
as an order of magnitude. For appropriate choices of the factorization scale, the difference 
can be reduced significantly, leading to more compatible results within the perturbative 
uncertainties, see the discussions in Refs. [21, 23-25, 28, 66]. As discussed already, the 5FS 
and 4FS possess different merits, and predictions that combine the advantages of both are 
highly desirable. The current combined values by the LHC Higgs Cross Section Working 
Group [67, 68] are obtained using the Santander matching prescription [30], which amounts 
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to a weighted average of the cross sections obtained in the two schemes. In contrast, our 
predictions here are derived from a consistent held-theory setup, and can thus be regarded 
as a definite improvement over the currently used prescription. 

We start in section 4.1 by extending the EFT result of section 2 to the case of two 
incoming protons. In section 4.2, we give details about the practical setup of our results. 
In section 4.3, we discuss our procedure to obtain robust estimates of the perturbative 
uncertainties from separate variations of the and iim matching scales. In particular, 
we discuss the profile scales and variations for the matching scale In section 4.4, 

we present our result for the bbH cross section as a function of the 5-quark mass. This 
serves as a validation of our matching procedure, confirming that our approach satisfies all 
the required properties. There, we also discuss the size of nonsingular power corrections 
suppressed by ra^jQ- Finally, in section 4.5, we present our final results at the physical 
6-quark mass for several Higgs masses and compare to the existing results obtained in the 
4FS, 5FS, and the Santander prescription. 


4.1 Extension of the EFT approach to hadron-hadron colliders 

The simplicity of the EFT framework presented in section 2 for DIS makes it possible to 
straightforwardly extend the setup to the case of two incoming protons. This is certainly 
not the case for any of the schemes discussed in section 3, whose consistent generalization 
to hadron-hadron collisions can be highly nontrivial. We first point out that the evolution 
of the quark operators and the matching at /im are identical and therefore the evolved 
PDFs of eq. (2.24) are the same. Of course, the matching at the hard scale is different. 
For m;, Q we have 


ObbH{mH,rnb) = Dij{mH,mb, 

while for rub <C Q we find® 


(4.1) 


— c c 

ObbH{mH,mb) = Cij{mH,'mb,HH)0\ 




[ 5 ], 




l[5]. 


+ Cbb{mH,^im) + oyMoyM 


[ 5 ]/ 


l[5]. 




(4.2) 


where we have introduced a bbH operator ObbH, and mn is the Higgs mass. Note that 
we have used the identities Cbk = C'bfc = C'fcb = Cj^i, and Cbb = Cn- These 

two results are the straightforward extensions of the results in eqs. (2.2) and (2.17), where, 
as discussed in section 2.4, we have made the choice to absorb power corrections into the 
coefficients for the light channels. 

®In this section we restore the difference between b and b, and we omit the convolution symbol (8) for 
ease of notation. 
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As in the DIS case, in our order counting we take the bottom PDF to be an object of 
order as- As discussed in section 2.3, this is appropriate for a hard scale of the order of 
the Higgs mass, (more generically for fin < ITeV). Therefore, up to NLO, the 

fixed-order result our cross section matches into for ij-h is given by (omitting the 

arguments for simplicity) 

LO(FO, 4F) a= 

NLO (FO, 4F) + ^^ 

+ ... (4.3) 

For Hm < fJ-H the resummed and matched cross section is written as^® 

LO+LL a= alciffJj + aH^cil'>fJg+ 2C^'> f\h 

NLO+NLL + fjj + hh + aH2C'£ hh 

+ ..., (4.4) 

where as in eq. (2.62) we have left implicit the strict expansion of the products of effective 
PDFs and coefficient functions. We notice that in both cases, using the perturbative order 
counting introduced in section 2.3, LO(+LL) is in fact order and NLO(+NLL) includes 
the order corrections. In section 4.4, we discuss the implementation and results of a 
strict expansion of eq. (4.4) as well as a more practical implementation keeping higher-order 
cross terms and using standard 5F PDFs. 

The 4FS result corresponds to the result of eq. (4.3) used for all scale hierarchies and 
where the decoupling scheme is used for the 6-quark renormalization of a*. The 5FS result 
on the other hand corresponds to the massless limit of eq. (4.4), replacing with 
and with the perturbative order counting performed only on the coefficient functions (i.e. 
assuming the bottom PDF of order 1). This has the expansion 

LO (5F) u = 

NLO (5F) + aj,2C(f/f/f + /f' 

NNLO (5F) + (2C^'> + 2Cj?) /f'/f ^ '/f + '/f' 

+ ..., (4.5) 

where CW are the massless coefficients (namely the massless limit of C^j )■ In figure 8 
we illustrate the different countings diagrammatically. This highlights that one can regard 
our results as a resummation-improved 4FS result. 

The massless coefficients required to reaching NLO+NLL 

accuracy in our result, eq. (4.4), are the same as those of the massless 5FS computation 
and can be found explicitly in ref. [28]. Trivially extending eq. (2.60) to the case of two 

^°For ease of notation, we do not distinguish between bottom and anti-bottom PDFs, and also on whether 
they come from one or the other proton, and compensate for this with numerical factors. 
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light-light 


light-heavy 


heavy-heavy 



LO-hLL 



NLO-hNLL 


5F NNLO 


Figure 8. Sample diagrams appearing in the computation of the Higgs production cross section in 
association with b quarks. Diagrams are grouped according to the different countings adopted in 
our resummed result, eq. (4.4), and in the 5FS result, eq. (4.5). The 4FS counting coincides with 
the resummed counting in the fixed-order limit where only the diagrams in the first column are 
considered. 


“ (2) ~(3) 

initial-state legs, the matching coefficients Ch' and Ch ' can be written as, 
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(4.6a) 

(4.6b) 


(4.6c) 

(4.6d) 

(4.6e) 


The MS massive coefficients can be obtained from the decoupling-scheme coefficients as 
described in Appendix A: 


^MS,(2) ^ ^(2) 






^MS,(3) ^ ^(3) _ 2^^ 








(4.7) 
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We have implemented analytic expressions for the coefficients D ^-in an in-house code and 
extract the numerical result for D\- from Madgraph5_aMC@NL0 [69], after generating the 
process pp —)• bbH at NLO. We have explicitly checked that our implementations, including 
pole scheme to MS scheme changes for in Dij, exactly reproduce the inclusive results of 
the bbhOnnlo code [28] and of the recent bbH studies of ref. [65]. 

4.2 Setup 

Here we summarize the set of input parameters we use to produce the results of sections 4.3 
and 4.4. Unless indicated otherwise we always use the setup detailed below. 

Collider energy We provide predictions for the LHC at ^/s = 8 TeV. 

PDFs We have created PDF sets using a modified version of APFEL [70] for the evolution 
from a fixed low scale where the parametrization of a known PDF set (MSTW2008) 
has been used. The main reasons for our modifications were the implementation of 
a general value for the threshold matching scale pm as well as the generation of the 
effective PDFs required in a strict expansion of eq. (4.4). Further details are 
given in Appendix B. 

Higgs mass We use niH = 125 GeV as default. 

Bottom mass For all results where the bottom mass is fixed to its physical value, we use a 
pole mass of rrn, = 4.75 GeV for the kinematic mass scale that enters in the 4F matrix 
elements and in the low-scale matching coefficients Mij. For the Yukawa coupling we 
use the MS mass mb{rnb) = 4.16 GeV as input, see also below. The use of different 
bottom masses for the Yukawa coupling and in the matrix elements is not unsual - 
this has been the setup of the 4FS calculations of Ref. [24, 25]. What is different 
to previously used setups (and to the LHGHXSWG) is that the two values we use 
are not related to each other via a one-loop conversion. This is not a problem, since 
the two perturbative series in which they enter are unrelated.What is relevant in 
our case is that we consistently use common values in both the resummation and 
fixed-order parts of the calculation. The numerical values above are chosen to have 
reasonable physical values and to enable an as consistent as possible comparison with 
the default 4FS and 5FS results. 

In our results where we vary mb to study the dependence on the bottom mass, mb 
and mb(jnb) are varied consistently, with the conversion between the two at one loop 
as required for our NLO calculation. 

Yukawa couplings All the Yukawa couplings are set to zero except the bottom quark 
Yukawa, Yb- The bottom Yukawa is renormalized in the MS scheme and its run¬ 
ning is set to 4 loops. In our numerical studies we always evaluate it at the hard 
scale ph, which is the appropriate scale for the resummation of large logarithms 

the future, a better approach would be to replace the pole mass in the threshold corrections by a 
proper short-distance mass scheme with a well-defined conversion from the MS scheme. 
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ln{fiH/mb) associated with the bbH vertex and hence leads to better perturbative 
convergence [71]. 

Bottom loops Contributions to Higgs production where the Higgs couples to a closed 
6-quark loop are usually included in the gluon-fusion cross section, since their most 
important effect is due to the interference of the bottom loop with the top loop. As 
usual, we exclude these contributions from our bbH computation, such that the result 
has no double counting with the gluon-fusion cross section. Our result still includes 
bottom loop contributions, but only in diagrams with two bottoms in the final state, 
(not included in the gluon-fusion cross section) as part of the NLO correction to the 
gg channel in our result. 

Yb ■ Yt interference In our results we neglect the interference contribution proportional 
to Yt ■ Yb by setting 1) = 0. In the SM, this correction is known to be important 
and reduces the inclusive 4FS NLO cross section by roughly 10% at the LHC for 
niH = 125 GeV [24, 25, 65], while in BSM scenarios with large tan/3 its relative 
contribution can be much smaller. This interference has been computed in the 4FS 
where it first enters at NLO via diagrams containing a top-quark loop, whilst in the 
5FS up to NNLO this interference does not contribute [28]. For comparisons between 
4FS and 5FS predictions it is often preferred that the interference terms are dropped 
[30, 72] since the latter are not present in the 5FS. To better compare with the 
results in the literature we also make this choice here. However, we emphasize that 
the Yt ■ Yb terms can be straightforwardly and consistently included as an additional 
nonsingular fixed-order piece in our result. To do so, we can simply allow for a 
nonzero top Yukawa in the fixed-order coefficients Dij. No changes to the resummed 
part of our result are required at the order we are working. 

4.3 Scale dependence and theory uncertainties 

In this subsection, we discuss in detail the perturbative uncertainties in our results. We 
begin by looking at the hard scale dependence. We fix mb to its physical value, set fim = 
mb, and plot in figure 9 the cross section obtained according to the 4FS (at LO and 
NLO, obtained using the code of [65]), the 5FS (at LO, NLO, and NNLO, obtained using 
bbhOnnlo [28]) and our result (at LO-I-LL and NLO-I-NLL). 

As expected, a clear reduction of the scale dependence is observed in all results when 
moving to higher orders. We also notice that the patterns of scale dependence of the 4FS 
(green dashed) and 5FS (blue dotted) results are opposite to each other with the former 
decreasing and the latter increasing with increasing fijj (except at NNLO). This is due to 
the fact that at LO the scale dependence is dominated by as for the 4FS result (which 
clearly increases at small scales), while for the 5FS result it is driven only by the bottom 
PDF, which vanishes at the bottom threshold and therefore drops rapidly as the scale 
decreases. Therefore, over a wide range of hard scales the two results differ significantly. 

In contrast, the framework we have presented in section 2 leads to cross sections (red 
solid) that are less sensitive to the choice of the hard scale, even at LO. The reason behind 
this is a large compensation between the contributions from the 66, bk, and ij channels. This 
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Figure 9. bbH cross section under hard-scale variation. 

is due to the fact that each 6-initiated contribution compensates the collinear subtraction in 
a gluon (or light quark) initiated contribution and close to the heavy-quark threshold these 
terms are all of the same order. This leads to a scale dependence in the (N)LO-|-(N)LL 
results that has a similar pattern to that of the unresummed 4FS result, however the 
resummation of collinear logarithms significantly stabilizes the dependence on fin- As with 
the 4FS, the 5FS results also have a greater dependence on fin compared to our resummed 
results. The reason for this is that the 5FS predictions adopt a standard perturbative 
counting and thus the compensation observed in the EFT results is not present. 

Additionally, figure 9 illustrates that a smaller scale fin ~ /4: leads to a more stable 

perturbative expansion for all the results, and also leads to better agreement between 
the different approaches. The reason for this has been studied in ref. [23] by a careful 
investigation of the actual size of the logarithms that arise in the 4FS prediction. 

Next, we discuss the choice of fim and its associated perturbative uncertainties. For 
this purpose, it is important to identify the kinematic region where the resummation is 
important and where it must be turned off. To this end, in the left plot of figure 10 we 
show the fixed NLO result and its decomposition into singular eq. (2.53) and nonsingular 
eq. (2.54) contributions, for a fixed value of the Higgs mass mn = 125 GeV and as a 
function of the bottom mass. In this plot we vary the bottom mass but have divided 
the cross sections by the bottom Yukawa coupling to better highlight the perturbative 
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Figure 10. Comparison of the magnitude of the singular, nonsingular, and full NLO cross sections 
when varying mb (left) and profile scale variations (right). See text for further details. 

structure. 

In the mb —t 0 region, the singular terms clearly dominate, while the nonsingular 
corrections are suppressed by at least an order of magnitude and tend to zero. This is the 
resummation region, where the canonical choice = ^6 is appropriate to resum the large 
logarithms ln(mfe/mjy) in the singular corrections. 

With increasing mb the singular contribution starts deviating from the full result, 
crosses it at around mb ~ 30 GeV ~ mjy/4, and becomes much larger than the full result 
in the large-rub region. This large-mb region corresponds to the fixed-order region, which 
exhibits a delicate balance between singular and nonsingular contributions, with a large 
cancellation between the two yielding the full result. This means that the distinction 
into singular and nonsingular is meaningless here. To not spoil this cancellation it is 
imperative that the resummation is switched off completely, which is done by taking = 
fiH- The fixed-order region starts at mb > mjy/d, where the magnitude of both singular and 
nonsingular is larger than the full result, so there is clearly an 0(1) cancellation between 
them. We have verified that this pattern holds at both LO and NLO and upon variation of 
the hard scale in the range mjy/16 < < mH- We can therefore safely take mb mn!^ 

as the point where we should turn off the resummation for any conhguration we might 
consider. 

A smooth transition between the canonical value = w-b in the resummation region 
and /im = in the fixed-order region is achieved by using prohle scales [56, 57], where 
the scale is promoted to a function of mb, which smoothly interpolates between these 
two limits. The use of profile scales is a common practice when performing resummation 
in EFTs based on RGEs. Following refs. [58, 59, 73], we choose different sets of profiles 
that allow us to separately estimate fixed-order and resummation uncertainties, which in 
the end are added in quadrature. 
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For our central scales we use 


fXH = m/f/ 4 , 
with the profile function 




f ( rrib \ 


X, 


if 0 < X < Xi 

X + 

(2-X2-X3){x-xi)^ 

if Xl < X < X2 

2(x2—Xi)(x3—Xi) ’ 
{2-Xi-X2){x-X3)^ 

I - 

if X2 < X < X3 

2{x3-xi){x3-X2) ’ 

1 , 


if X3 < X 


(4.8) 


(4.9) 


and we have chosen the appropriate values of {xi = 0.3, X 2 = 0.65, X 3 = 1.0}. In this 
way, the resummation slowly turns off as increases, becoming completely switched off 
for mb > m/f/4, which corresponds to the point identified above. The right-hand plot of 
figure 10 illustrates these profile functions: the solid green curves correspond to eq. (4.8) 
as a function of mb for the hard scale choices hh = {mn/8,mH/4:,mH/2}. Note that 
at small mb the standard scale Hm = mb is recovered for the central profile scale with 
Hh = m///4. The hard scale variation by a factor of two leaves the ratio fim/fJ-H fixed and 
therefore does not change the resummation. At the same time for large mb it recovers the 
usual fixed-order scale variation. Hence, we use these variations to estimate the fixed-order 
uncertainty Apo- 

The variation of the central prohle, while keeping the hard scale fin fixed, is performed 
by multiplying the central profile by a factor. 


f4^''^{mb,HH,mH,a) = /vary ^mimb, , mn) 

_ fa ( \ j, ( mb \ 

I^H /vary J Vm^/d J ’ 


(4.10) 


where a e [— 1 , 1 ] and 


/vary ( 3 :) — < 


2 ( 1 -!^,, 


1 + 2 ( 1 - 

1 , 


X3 


if 0 < X < f 
if ^ < X < X 3 
if X3 < X. 


(4.11) 


The multiplicative factor /vary( 3 ^) tends to 2 in the limit x —)• 0 and tends to 1 in the limit 
X —)> X 3 (as before, we use X 3 = 1). The effect of this factor (when varying a G [—1,1]) 
is to vary the arguments of the resummed logarithms in the small mb region by a factor 
of two, while keeping the hard scale fixed. Hence, we can use these variations to estimate 
the resummation uncertainty Aresum- In the limit x —)■ X 3 (or mb —)■ fJ-n) the effect of 
this variation tends to zero, as it must, and thus the resummation uncertainty vanishes in 
the fixed-order result as it should. In the transition region between the resummation and 
fixed-order regions, this variation effectively captures the uncertainty in the transition. In 
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the right-hand plot of figure 10 the yellow band enclosed by the dotted green curves shows 
the effect of this variation on the central profile fin = 

A key advantage of the setup we have discussed above is that it provides a concrete 
way by which to estimate the theoretical uncertainties. Our result for the cross section is 
obtained via a two-step matching procedure and the variation of the scales at which this 
matching is performed is a natural way to arrive at a realistic error estimate. To obtain 
our estimate of the total theoretical uncertainty, we take Apo and Aresum as the maximum 
variation among each of their respective profile variations. We then obtain Atot by adding 
the two in quadrature, 


^tot — ^FO + ^resum- (4-12) 

We emphasize that since all variations we perform amount to variations of scales (albeit 
more intricate than standard scale variations), the resulting perturbative uncertainties 
ApO) Aresum and Atot decrease when increasing the perturbative order of a calculation. 

Finally, we note that profiling the scale Hm is nontrivial for general values of mb and 
^H- Since fim corresponds to the scale at which PDFs are matched from a theory involving 
bottom quarks to a theory with no bottom quarks, each point of a profile function for 
corresponds to a different PDF set (with as the bottom threshold). To produce the 
rufj-variation plots in section 4.4 we have produced 20 PDF sets for each of the five profiles 
in figure 10. For the results at the physical value presented in section 4.5, we are always 
in the canonical region [x < xi in the profile function), which means we are only required 
to generate PDF sets for the values £ {0.5mft, m;,, 2mfe}. 

4.4 Cross section and power corrections as a function of mb 

In this subsection we study the cross section as a function of mb- The reason for this is 
to confirm that the result obtained in the framework presented in this paper does indeed 
smoothly interpolate between resummation and fixed-order regions. It also serves as an 
important validation of the method we employ to estimate uncertainties. In the left-hand 
plot of figure 11 we show the LO-I-LL (dashed dark green) and NLO-I-NLL (solid navy) 
cross sections including error bands (green and blue bands respectively). We also plot 
central values for the associated fixed-order cross sections at LO (dotted dark green) and 
NLO (dashed navy). The right-hand plot of figure 11 displays the relative size of the total 
LO-I-LL uncertainty (green band) and of the NLO-I-NLL resummation (light blue band) 
and total (navy band) uncertainties. We emphasise that the results in figure 11 are from 
an implementation of the strict expansion of the cross section in eq. (4.4). 

The first feature to point out is that at large mb both LO-I-LL and NLO-I-NLL results 
tend to their fixed-order counterparts (i.e., tend to the LO and NLO cross sections). This 
clearly shows that the framework we have introduced indeed fulfills this desired property 
in the limit of large mb- The fact that this transition occurs smoothly is a natural result 
of the strict expansion used here (that is, there are no higher-order cross terms in our 
result that might spoil the full cancellation between resummation pieces). The smooth 
transition between the low and high mb regions is also a direct consequence of the order 
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Figure 11. (N)L0+(N)LL cross section for bbh as a function of nib (left) and relative uncertainties 
(right). A strict expansion of the cross section eq. (4.4) has been performed. 


counting we adopt. If we were not to regard the 6-quark PDF as an order as object, then 
the strict expansion we have performed would not be possible and a discontinuity would 
be present in the region of mb ~ fin arising from the higher-order terms. Finally, the 
smooth transition to the fixed-order results indicate that our method for including the 
power-suppressed 0{ml/m‘jj) terms to the strict EFT result works perfectly, and that as 
we have discussed in section 2.4 it is indeed the case that (at least to the order we work to) 
it is possible to consistently include all power-corrections present in the fixed-order result. 

Regarding the estimates of the perturbative uncertainty, figure 11 reveals that the 
error bands we assign are indeed reasonable and robust over the full range of mb, with the 
NLO-I-NLL band fully contained within that of the LO-I-LL. The right-hand plot of figure 11 
indicates that the total uncertainty is dominated by the fixed-order scale uncertainty in 
the large-rub limit, and the resummation uncertainty vanishes, as it should, in the limit 
fJ-m However, with decreasing mb we see that the resummation uncertainty becomes 

nonnegligible, forming an important component of the total error. The total uncertainty 
is of the order of 12-14% over most of the range of mb considered here, and grows as mb 
is increased beyond the scale fin- 

Our result also allows us to consistently quantify the size of power corrections of 
0{m'l/m?^). This relies on the observation that in the small mb limit, due to the vanishing 
of the nonsingular contributions, our result essentially becomes a re-arranged 5F computa¬ 
tion. This means that all terms required to obtain a consistently matched prediction can in 
principle be extracted from a calculation that sets mb = 0 from the outset. The comparison 
between such a re-arranged 5F calculation and the result where the power corrections are 
included allows us to study the size of the latter. To illustrate this, in figure 12 we compare 
the LO-I-LL prediction where power corrections in mb have been included (solid blue) to 
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Figure 12. A comparison of the LL+LO cross section where nonsingular corrections ^ 
are included (solid dark green) and have been set to zero (dotted gray). The latter cross section 
is essentially a re-arranged 5F computation. The lower panel indicates the size of such power 
corrections through the ratio of the two cross sections. For reference we have included the LO 
(dotted blue) and NLO (dashed blue) 5F predictions in the upper panel. 


the same prediction made with strictly massless coefficient functions (dotted gray).^^ It is 
clear that in the small limit the nonsingular terms are unimportant and 

that the matched result can simply be constructed, with negligible errors due to missing 
power corrections, from massless coefficient functions. This argument indicates that should 
S-ACOT or FONLL with standard perturbative counting be applied to the case of bbH, 
then the resulting cross section will likely be almost the same as that of the 5F prediction. 

It is also apparent that these power-corrections do increase in importance, their size 
exceeding 10%, for mb > 10 GeV. Therefore, in such parameter regions including them is 
vital for a faithful description of the cross section. In figure 12 we have also plotted the 5F 
LO and NLO results, which deviate visibly from the LO-I-LL result as mb grows, indicating 
that in such regions a massless 5F prediction becomes an inadequate description of the 
process. 

Finally, we make some brief comments regarding the difference between the cross sec¬ 
tion obtained under a strict perturbative expansion of eq. (4.4) (or equivalently eq. (2.61)) 
compared to that obtained by not expanding the two sets of matching coefficients. As 

^^Given that the term is not required to make a prediction at NNLO for bbH in the 5F scheme, 
this ingredient in the massless limit is not known analytically and therefore we could not make the same 
comparison for the NLO-I-NLL result. Nevertheless, by construction, exactly the same pattern that is 
observed for the LO-I-LL result is expected to hold for the NLO-I-NLL result. 
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Figure 13. NLO+NLL cross section for bbH using the nonexpanded implementation of eq. (4.4). 
The solid and dotted curves correspond to the expanded and nonexpanded predictions re¬ 
spectively and magenta, blue and green colours correspond to the hard scale choices = 

mn/S, mn/'i, 


mentioned earlier, not performing a strict expansion will generically lead to a discontinuity 
in the limit /im due to the non-cancellation of spurious higher-order interference 

terms. This is illustrated in figure 13 where we have plotted the NLO-I-NLL cross section 
predictions for three hard scale choices under a strict perturbative expansion (solid) and 
with no strict expansion (dotted), namely using standard PDFs ff'^. At large nib the solid 
and dotted curves display significant differences and in particular the latter showing dis¬ 
continuities for —)• In the small mb limit however, it is clear that the differences 

between expanded and nonexpanded approaches become much smaller. In particular, this 
means that the total uncertainty band in the region of physical 6-quark masses is basi¬ 
cally the same in the two approaches. This property can be used to greatly simplify the 
practical implementation in this region and we have exploited this to produce the results 
in section 4.5. However, it is also important to point out that in general it is only an ex¬ 
panded result, akin to that of eq. (2.49), that guarantees a consistent and smooth matching 
between resummation and fixed-order regions. These observations may well be important 
when considering heavy-quark initiated processes where the value of m/Q is not as small 
as in the setup we study here. 

4.5 LHC phenomenology 

Here we return to the phenomenologically relevant case of a physical bottom mass and 
consider the cross section as a function of the Higgs mass. The previous subsection con- 
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firmed that the framework we use leads to a cross section that consistently describes both 
resummation and fixed-order regions. It is of course of interest to compare in a meaningful 
manner our (N)LO-|-(N)LL predictions to other predictions available, namely predictions in 
the 4F and 5F schemes, as well as to the Santander matching prescription used to combine 
these two. The latter is a practical formula that combines the 4FS and 5FS predictions for 
the total inclusive cross section through a weighted average of the two [30] , 


^matched 


^4FS^^^5FS 
1 + U1 


where 



(4.13) 


This construction is such that the combined result tends to that of the 4FS when the 
collinear logarithms, ln(m///mb) are small, and to that of the 5FS when the logarithms 
are large (i.e., in the limit mh/mH —)• 0). The choice of the weight u is motivated by the 
fact that this leads to roughly equal weights being assigned for the 4FS and 5FS numbers 
around mu ~ 100 GeV, which is the region of ‘best’ agreement between the 4FS and 5FS 
predictions (see for example ref. [66] ). Beyond this motivation the choice of u is arbitrary 
and there is no strong theoretical argument preventing the choice of alternative weights 
or different ways of averaging the two cross section predictions. Moreover, the practical 
formula combines two predictions made using different PDF and rribijnh) inputs, which is 
somewhat inconsistent. The estimate of the uncertainty on a Santander matched prediction 
is given by the error band obtained by applying the formula eq. (4.13) to the upper and 
lower uncertainty curves of the 4F and 5F predictions. 

Regarding the set of inputs we use, we have chosen to stick as closely as possible to 
those used in the LHCHXSWG [67, 68, 72] and also those used in recent studies of hbH 
production [65]. Explicitly, we use the default MSTW2008 PDF sets, at the appropriate 
order for the LO, NLO and NNLO 5FS predictions, whilst we use the hxed-flavour nj = 4 
set for the 4FS predictions. For the 4FS and Santander matched results we explore the effect 
of using mb{fnb) / 4:. 16 GeV (i.e., a different rflbimb) than that used in 5FS predictions), 
as done by the LHGHXSWG.^^ The central choice of hard scale is = {mu + 2mb)/4, 
with mb = 4.75 GeV, and we vary this hard scale by a factor of two to obtain the fixed- 
order uncertainty. The errorbars for the 4FS and 5FS predictions are obtained by setting 
fJ-F = fJ-R = that is we do not consider fip 7 ^ fJ-R variations here. The bands for the 
Santander matched cross sections are obtained with the Santander prescription. 

In hgure 14 we plot the 5FS (blue points) and 4FS (green points) cross sections, the 
(N)LO-|-(N)LL matched predictions (red points) as well as the Santander matched cross 
sections (brown and magenta points) for mn G {110,125,140} GeV. The error bands for 
the (N)LO-|-(N)LL predictions have been obtained as discussed in eq. (4.12). We note 
here that since the ratio mfe/(m/f/4) < 0.3, the profile region we are in is actually always 
linear, namely Umimt, rnn, fJ-n) oc mb- Gompared with the green 4FS NLO point, the light 
green 4FS NLO point has been obtained by setting fnb{fnb) = 4.34 GeV. This (somewhat 
artificially) shifts the cross section upwards by ~ 10% by increasing Yb{fiH)- The magenta 

^®The reason the choice mb{fnb) = 4.34 GeV 7 ^ 4.16 GeV is made in some 4FS predictions is that this 
MS mass corresponds to a 1-loop conversion of the used pole mass. However, by using a fixed pole mass as 
input one reintroduces the pole mass renormalon ambiguity into the cross section through Vj,. 
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Figure 14. A comparison between all available cross section results for bbH. 


Santander point has been computed with the 4FS component using this increased value of 
bnbijnb) and is therefore seen to be higher than the point that uses a consistent MS mass 
in both 4FS and 5FS components. We have checked that the results of figure 14 are fully 
consistent with those in the literature. 

As seen in the previous subsection the error band of the NLO+NLL predictions is 
contained within that of the LO+LL. The LO+LL uncertainty band is quite wide (and 
larger than the 4FS or 5FS LO bands) — something that is to be expected from a LO 
prediction and a maximal variation of the two matching scales, eq. (4.12). Rather than 
being of concern, this observation gives us confidence that we do not underestimate the 
inherent uncertainties present and furthermore allows us to trust the size of the NLO+NLL 
band. We additionally observe that most of the LO+LL scale uncertainty is driven by the 
variation, while at NLO+NLL the uncertainty band is dominated by the variation. 
On the other hand, the 4FS displays a less significant reduction in its uncertainty band at 
NLO, and with its large correction shows relatively poor perturbative convergence due to 
the presence of unresummed logarithms. The 5FS band shrinks more visibly with increasing 
order, however the NNLO central value lies outside the NLO band.^'^ 

note that the 5FS bands (and to a much lesser extend the 4FS bands) would be larger if we had 
considered fj,R ^ fj,F, potentially improving the convergence pattern described above. We have not done 
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The (N)LO+(N)LL results lie significantly higher than their respective 4FS (N)LO 
counterparts. This is due to the resummation contained in the former results. We also 
notice that the NLO+NLL results lie slightly above the NNLO 5FS predictions. The reason 
for this is that the former results contain the (positive) effects of light channels {gg, qg 
and qq) at whilst the latter results contain (negative) two-loop corrections to the bb 

channel (see figure 8). 

By construction the Santander matched result lies between the 4FS NLO and 5FS 
NNLO predictions (irrespective of the precise inputs for mb (mb) used in the 4FS calcula¬ 
tion). Our NLO-I-NLL result is therefore higher than the Santander matched results, and 
specifically we find an increase of 6% with respect to the magenta LHCHXSWG Santander 
point, and of 12% with respect to the brown Santander point (which is directly comparable 
to the NLO-I-NLL result). We also notice that the error bands of the NLO-I-NLL prediction 
are roughly of the same size as those obtained through Santander matching indicating that 
the size of the latter is likely to be realistic. There is a sizeable overlap in the uncertainty 
bands of the Santander and the NLO-I-NLL predictions, however the central Santander 
points lie towards the bottom edge (magenta) and outside (brown) of the NLO-I-NLL error 
band. We also point out that while our NLO-I-NLL result is stable upon hard scale varia¬ 
tion, the Santander matched result would be significantly smaller for larger hard scales, as 
is clear by inspection of figure 9. 

In figure 15 we present the same comparison as in figure 14 for larger Higgs masses, 
niH = 300 GeV and mn = 500 GeV. The overall pattern observed for a light Higgs 
remains qualitatively unchanged. The NLO corrections in the 4FS have grown noticeably, 
such that the 4FS NLO results are now outside the 4FS LO bands, which is likely due 
to the unresummed logarithms, which get larger at higher Higgs masses. While at lower 
niH the central values of the 4FS NLO and 5FS NNLO are relatively close to each other 
and have overlapping uncertainty bands, at higher mn they are further apart and have 
only barely overlapping uncertainty bands. Gonsequently, the Santander average becomes 
even less reliable. Encouragingly, the LO-I-LL and NLO-I-NLL results at large niH continue 
to display the good perturbative behaviour and convergence pattern present at lower mn 
values. They also remain systematically higher than the Santander matched predictions. 

Finally, we briefly comment on the effect of including the known term^^ (i-e., the 
two-loop corrections to the 66-channel), which is formally of higher order in our approach. 
This is the only known contribution that we have not included in the NLO-I-NLL result. 
The NLO-I-NLL result for mn = 125 GeV shown in figure 14 is unlo+nll = 0.224 ± 
0.021 pb. With the addition of the higher-order C^T-term this cross section becomes 

1 ^( 2 ) = 0.211 ± 0.010 pb, to be compared with the NNLO 5FS cross section 
00 

of o' 5 f,nnlo = 0.209 ± 0.010 pb. Glearly the addition of this higher-order term reduces 
the NLO-I-NLL cross section and additionally reduces its uncertainty band, bringing both 
central value as well as the size of the error bands closer to those of the 5FS NNLO 


this in order to directly compare to our method of estimating uncertainties. 

Additionally, we include the lowest order 66-channel term appearing at the same order, which 

however gives a negligible effect. 
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Figure 15. Same as Fig. 14 but for large Higgs masses. 

result. This term is likely to be important at very high scales ^ 1 TeV, where the 5FS 
perturbative counting is expected to be more appropriate. However, we feel that including 
this term for a SM Higgs is rather ad-hoc given that there are a number of additional 
terms that would also contribute at the same order, but which are not known and have 
not been included here. Furthermore, the sizeable reduction of the error bars is slightly 
discomforting given that we have no control on the effects of these missing terms. We note 
that the error bars of the original NLO-I-NLL cross section nicely cover the effect of adding 
the term, which provides further support that the uncertainty bands presented are 
indeed reasonable. 

5 Conclusions 

We have presented a systematic EFT setup to derive heavy-quark initiated cross sections 
at hadron colliders. Our framework includes the resummation of potentially large loga¬ 
rithms ~ ln{mh/Q). Furthermore, it consistently includes power corrections ~ 
reproducing the full fixed-order (4FS) result. As such our final result gives predictions that 
are accurate in both of the limits mb Q and mb ~ Q as well as in the transition region 
in between. 
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Our result is obtained via a two-step matching procedure, and variation of the scales 
at which these two matchings are performed allows us to obtain a robust estimate of the 
perturbative uncertainties. The construction of the coefficient functions of our result bears 
several similarities with existing VFNSs for DIS. A key difference in our approach is the 
different perturbative order counting. In particular, we argue that it is more appropriate 
to count the effective 6-quark PDF, which is generated perturbatively at the scale fj-m, 
as a perturbative object of 0{as). This organization of the perturbative series leads to 
perturbatively stable results and allows for a smooth transition between fixed-order and 
resummation regions. The simplicity of the EFT approach for DIS makes its generalization 
to hadron-hadron collisions straightforward. 

We have applied our framework to the case of the bbH cross section. We first studied 
the cross section as a function of mf,, which served to demonstrate that our resummed and 
matched result satisfies all required properties. We then presented numerical results of 
phenomenological interest for the LHC. We have compared our results to the 4FS and 5FS 
result, as well as the Santander prescription, which combines the two results by taking a 
weighted average. Since our predictions are derived from a field-theory setup, consistently 
combining the 4FS and 5FS limits, it can be regarded as a definite improvement over this 
prescription. At NLO-I-NLL, we find a slightly reduced perturbative uncertainty compared 
to the Santander average. The Santander central value is lower than the NLO-I-NLL result 
by about 12% for mn = 125 GeV and lies outside our uncertainty band when consistent 
MS masses are used in both the 5FS and 4FS ingredients of the Santander result. The 
difference is reduced to 6% when using a larger MS mass in the 4FS result, as used by the 
LHCHXSWG, with the central value lying at the lower edge of our uncertainty band. We 
observe that the NLO-I-NLL result is stable upon variation of the hard scale, unlike the 
Santander matched results which would vary significantly under such variation. 

The framework presented in this paper can be straightforwardly applied to other pro¬ 
cesses involving initial-state 6-quarks, for example single-top or E -|- 6-jet production. Ad¬ 
ditionally, extending the framework to study more differential observables of interest, such 
as the transverse momentum of the Higgs in 66iL-production or jet-vetoed cross sections is 
possible. In making the calculation less inclusive, more scales appear in the perturbative 
expansion, which can be efficiently dealt with in an EFT setup. Finally, it would be very 
interesting to adopt our perturbative counting in DIS in the context of PDF fits, where 
the used variable flavor number scheme typically plays a central role in determining the 
accuracy and the goodness of the fit itself. 


Note added: While this paper was being finalized ref. [74] appeared, which obtains the 
bbH cross section in the FONLL approach. As discussed in section 3, the FONLL approach 
adopts a different perturbative counting to the one we have presented here. In particular, 
the result of ref. [74] is computed at FONLL-A accuracy, that is, it combines the 4FS LO 
result with the 5FS NNLO result. In comparison to our NLO-I-NLL result, it does not 

include the 4FS NLO contributions from Cgg , , and Cqg . However, it does include the 

(‘ 2 ') ( 2 ) 

and 5FS NNLO terms, which are higher-order terms in our approach (and whose 


- 44 - 


impact on our result is discussed at the end of section 4.5). As a result, the FONLL-A 
result of ref. [74] is very close to the 5FS NNLO, as one might expect, since the difference 
to the latter is the inclusion of the numerically small 0 (ag power corrections from 

the 4FS LO result. Finally, ref. [74] does not include an estimate of the resummation and 
matching uncertainty (analogous to our /j,m variation). 
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A Renormalization 

In this appendix we briefly discuss aspects of renormalization important to our work. In 
the MS scheme ultraviolet (UV) divergences associated with quark lines are subtracted in 
the same way for all quarks, independently of their mass (mass doesn’t play a role in the 
UV). However, a direct consequence of the MS scheme is that a* runs with n/ = 6 flavors, 
irrespectively of the energy scale. A more physical renormalization scheme for heavy quarks 
is the so called CWZ or decoupling scheme [75], where all “light” quarks are renormalised 
with MS counter terms, while UV divergences associated with the heavy-quark loops are 
subtracted at zero momentum. This ensures decoupling of the heavy quarks at energies 
much smaller than their mass. Since the concept of light and heavy depends on the actual 
scale, in practice a variable flavor number renormalization scheme is used, where the number 
of active (light) quarks renormalised in MS depends on the hard scale and changes at the 
crossing of the heavy quark thresholds. We have used this in Sect. 2 when treating the 
fixed-order and resummation regions differently regarding the renormalization of 6-quark 
loops. 

At the threshold scale p, = prm matching conditions relate the value of as above and 
below that scale. Denoting with a superscript the number of flavors used in the evolution 
of Us, we have 

afHhm) ^ Tf ^ 

^ - 3 

where m is the heavy quark pole mass. These conditions are currently known through 4 
loops [76], but for our applications we just need the 2 loop expression of eq. (A.l). 

A generic observable F{as) can be written as a perturbative expansion equivalently in 
both schemes, 

F{as) = j;aW\/r)FW(^)(/r) = Y, , (A.2) 

k k 
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and to all orders they are identical. The relation between the coefficients whose 

heavy quark UV divergences are renormalised in the decoupling scheme, and renor¬ 

malised in MS, can be simply obtained by using eq. (A.l) to write in terms of in 
the first sum (or viceversa), re-expanding and matching order by order. 

B Construction of PDFs with variable threshold 

In this work we have considered PDFs in which the heavy quark thresholds are not hxed 
to the heavy quark mass, but can vary. Moreover, the effective PDFs defined in eq. (2.48) 
are required with mixed evolution and matching accuracies. If we consider the “universal” 
PDFs as those at a small scale //Aj then the dependence on the threshold and the details 
of the order at which each ingredient is retained are all in the perturbative evolution. 
Typically PDFs are used through the LHAPDF library, where evolved PDFs at any (available) 
scale are built from interpolation grids, previously created assuming a particular evolution 
with specific heavy quark thresholds. 

For our purposes, performing the evolution each time picking the desired value of the 
heavy quark thresholds seems advisable. However, this approach faces speed problems, 
since performing the evolution is much more time consuming than interpolating a grid. 
Therefore, for the work in this paper we have created LHAPDF grids with different choices 
for the 6-quark threshold, and with the required combinations of the perturbative 
ingredients at different orders. To do so, we have used the public code APFEL, which 
has the ability of creating grids after performing its own evolution. To accommodate the 
possibility of choosing a threshold different from the heavy quark mass, we have modihed 
the code adding the //m-dependent matching conditions through NNLO from Ref. [53]. 
Furthermore, we modihed the code to produce PDF grids with the required combinations 
of orders of the matching coefficients, Mij. Practically, we proceeded as follows: 

• We start with a central member of a public PDF set, namely MSTW2008 with 
as{mz) = 0.1171, both at NLO and at NNLO. The value of the bottom pole mass 
used in this work has been taken to be nib = 4.75 GeV for consistency with the 
chosen PDF set. 

• We compute all the required PDFs as well as a* from this set at an initial scale 

fJ-o = = 1.4 GeV. 

• We use APFEL to perform the forward evolution and create the corresponding grids 
for the different setups we are interested in. 

Note that we choose to perform the evolution at (N)NLL, starting from a (N)NLO set, for 
producing the PDFs that we used in our (N)LO-|-(N)LL results. While this is somewhat 
in constrast with the discussion in section 2.3, using the evolution at one higher order has 
two advantages. The first is that the order of the evolution and the highest order in the 
matching functions are consistent, so that when we do not use the strict expansion we can 
use standard PDFs, as explained below eq. (2.50). The second is that our final NLO-I-NLL 
result in section 4 is more directly comparable to the 5FS NNLO result. 
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Figure 16. Bottom PDF as a function of the scale at fixed x = 0.005 for different values of the 
bottom threshold at NLO (left) and NNLO (right). 

As an illustration of what the modified code is capable of, in figure 16 we show the 

f5l 

standard 5F bottom PDF for a fixed value of x = 0.005. We plot this as a function of 
the scale /r close to the bottom threshold for three different values of the threshold itself, 
Hm = rnb/2,mb,2mf„ as used in the phenomenology section 4.5. Changing the threshold 
of course shifts the value of the PDF at fj-m, and in particular makes this value nonzero 
at NLO (at NNLO it is already nonzero even for = TUb). The initial condition, shown 
as the gray dotted or dashed lines, is given by the product of 4F PDFs with the matching 
conditions. We observe that at NNLO the PDFs obtained with different valnes of the 
threshold are almost identical at a large scale, indicating that a NNLO cross section is 
only likely to have a mild dependence on fj,m- We also note that the initial condition itself 
behaves in a nice pertnrbative way at large scales, where we expect higher order corrections 
to be small, while it deviates strongly at smaller scales, where Og is larger and higher-order 
corrections are not negligible. 

C Massive kinematics in hadron-hadron collisions 

In the case of an incoming massless parton its fonr-momentum is considered to be a fraction 
of the four-momentum of the proton. In case of massive incoming parton this formulation 
wonld lead to a mass that scales with the momentum fraction, which is clearly inconsistent. 
To overcome this, the proper approach [5] is to nse light-cone coordinates, where momenta 
can be written as 

p= p^ = {p^ ±p^)/V2, (C.l) 

where p* are the usual Minkowski components. We choose to orient the beam axis along 
the third spatial direction. A collinear particle with mass m has pt = 0, and the mass 
can be expressed as = 2p'^p~. We can write the momenta of massless protons^® in the 

^®Tlie proton mass can always be neglected compared to its momentum at the LHC. 
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hadronic center-of-mass frame as 


Pi = (p+,o,o), P 2 = (o,p^,6) 


(C.2) 


with P^ = P 2 = sjS/2, and S = (Pi + P 2 )^ is the total invariant mass squared of the 
colliding protons. The collinear component of a parton’s momentum scales with the largest 
light-cone component of the proton, so two incoming partons with momentum fractions xi 
and X 2 have momenta 


Pi 



9 

2xiP^ 



P2 


ml 

2X2 P 2 ” ’ 


X2P2 ,0 


(C.3) 


The accessible values of xi ^2 are determined by kinematic constraints. Imposing the con¬ 
dition that the parton energy cannot be greater than the proton energy, we find the con¬ 
straints 






i = 1,2. 


(C.4) 


Since the masses 2 are negligible with respect to 5, these conditions reproduce the 
massless limit constraint 0 < xi ^2 < 1- The partonic invariant mass squared is given by 


/ n 2 r, 2 2 mlml 

s = {pi +P 2 ) = X 1 X 2 S + mi + ml + ^ ^ g - 

In order to produce a final state of invariant mass M the inequality 

s> 


(C.5) 


(C.6) 


must be satisfied. If we consider at least one of the two partons to be massless (say m 2 = 0), 
this inequality has a single solution X 1 X 2 > /S. This is already problematic, since 

for small invariant masses M —?• mi very small values of xi ^2 are accessible. The situation 
becomes worse if both partons are massive, as in the case of the subprocess hh ^ H. In 
this case there are two solutions of eq. (C.6), namely (setting mi = m 2 = m for simplicity) 


X 1 X 2 > 


m2 

45 



and 


X1X2 < 


m2 

45 



(C.7) 


where the second solution represents a new region of very small Xi that is inaccessible in 
the massless case. The physical interpretation is as follows. As the momentum fraction 
of a parton is reduced, its energy reaches a minimum at x* = mj/\/5 (where it is not 
moving in the center-of-mass frame), and then starts increasing upon further reduction of 
Xi- Therefore, at very small Xj, a parton is very energetic again thus making it possible 
to produce a high invariant mass final state. In this configuration both heavy quarks have 
become “anticollinear” with respect to their respective protons. 

It is clear from this simple kinematical argument that a massive extension of the stan¬ 
dard factorisation theorem cannot just work in its usual form in presence of two incoming 
hadrons (otherwise, there would be a huge contribution from unconstrained small-x PDFs). 
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When there is just one proton, as in DIS, the problematic configuration described above 
never takes place, the effect of the parton’s mass being a further restriction to the acces¬ 
sible values of x with respect to the massless case, and standard collinear factorisation 
works even in presence of massive partons [5]. In the hadron-hadron collider case, only 
a systematic expansion in the heavy quark mass such as the one presented in this work 
allows the description of heavy quarks in the initial state. 
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